diff --git a/CMakeLists.txt b/CMakeLists.txt index abf3ee350..182ad0939 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -100,8 +100,17 @@ target_sources( src/ASM/BDFMats.C src/ASM/BlockElmMats.C src/ASM/DomainDecomposition.C - src/ASM/dslbln.f - src/ASM/dsolv3.f + src/ASM/Ftn/dcros3.f + src/ASM/Ftn/dmuax3.f + src/ASM/Ftn/dnorm3.f + src/ASM/Ftn/dprint.f + src/ASM/Ftn/dslbln.f + src/ASM/Ftn/dsolv3.f + src/ASM/Ftn/errmsg.f + src/ASM/Ftn/q4_inv.f + src/ASM/Ftn/shape2.f + src/ASM/Ftn/shltri.f + src/ASM/Ftn/t3_inv.f src/ASM/DualField.C src/ASM/ElmMats.C src/ASM/Field.C diff --git a/src/ASM/ASMLagBase.C b/src/ASM/ASMLagBase.C index f38ef3c46..0c8264352 100644 --- a/src/ASM/ASMLagBase.C +++ b/src/ASM/ASMLagBase.C @@ -12,6 +12,7 @@ //============================================================================== #include "ASMLagBase.h" +#include "Vec3Oper.h" ASMLagBase::ASMLagBase (const ASMLagBase& patch, bool cpCoord) : coord(myCoord) @@ -20,7 +21,8 @@ ASMLagBase::ASMLagBase (const ASMLagBase& patch, bool cpCoord) : coord(myCoord) } -bool ASMLagBase::nodalField (Matrix& field, const Vector& sol, size_t nno) const +bool ASMLagBase::nodalField (Matrix& field, + const RealArray& sol, size_t nno) const { size_t nPnts = coord.size(); size_t nComp = sol.size() / nno; @@ -28,7 +30,7 @@ bool ASMLagBase::nodalField (Matrix& field, const Vector& sol, size_t nno) const return false; field.resize(nComp,nPnts); - const double* u = sol.ptr(); + const double* u = sol.data(); for (size_t iPt = 1; iPt <= nPnts; iPt++, u += nComp) field.fillColumn(iPt,u); @@ -48,7 +50,22 @@ Vec3 ASMLagBase::getGeometricCenter (const std::vector& MNPC) const } -bool ASMLagBase::updateCoords (const Vector& displ, unsigned char nsd) +double ASMLagBase::getBoundingBox (const std::vector& MNPC, + Vec3Pair& bbox) const +{ + bbox.first = bbox.second = coord[MNPC.front()]; + for (size_t n = 1; n < MNPC.size(); n++) + for (int i = 0; i < 3; i++) + if (double x = coord[MNPC[n]][i]; x < bbox.first[i]) + bbox.first[i] = x; + else if (x > bbox.second[i]) + bbox.second[i] = x; + + return (bbox.second - bbox.first).length(); +} + + +bool ASMLagBase::updateCoords (const RealArray& displ, unsigned char nsd) { if (displ.size() != nsd*coord.size()) { @@ -58,7 +75,7 @@ bool ASMLagBase::updateCoords (const Vector& displ, unsigned char nsd) return false; } - const double* dpt = displ.ptr(); + const double* dpt = displ.data(); for (Vec3& XYZ : myCoord) { XYZ += RealArray(dpt,dpt+nsd); diff --git a/src/ASM/ASMLagBase.h b/src/ASM/ASMLagBase.h index dc991cb09..db0424113 100644 --- a/src/ASM/ASMLagBase.h +++ b/src/ASM/ASMLagBase.h @@ -34,15 +34,17 @@ class ASMLagBase virtual ~ASMLagBase() {} //! \brief Direct nodal evaluation of a solution field. - bool nodalField(Matrix& field, const Vector& sol, size_t nno) const; + bool nodalField(Matrix& field, const RealArray& sol, size_t nno) const; //! \brief Returns the geometric center of an element. Vec3 getGeometricCenter(const std::vector& MNPC) const; + //! \brief Returns the bounding box and diameter of an element. + double getBoundingBox(const std::vector& MNPC, Vec3Pair& bbox) const; //! \brief Updates the nodal coordinates for this patch. //! \param[in] displ Incremental displacements to update the coordinates with //! \param[in] nsd Number of space dimensions - bool updateCoords(const Vector& displ, unsigned char nsd); + bool updateCoords(const RealArray& displ, unsigned char nsd); public: //! \brief Updates patch origin by adding a constant to all nodes. diff --git a/src/ASM/ASMbase.C b/src/ASM/ASMbase.C index 369f130b6..225daf07f 100644 --- a/src/ASM/ASMbase.C +++ b/src/ASM/ASMbase.C @@ -1622,6 +1622,19 @@ bool ASMbase::injectNodeVec (const RealArray& nodeVec, RealArray& globRes, } +bool ASMbase::findPoints (std::vector& points, + std::vector& locs) const +{ + size_t nFound = 0; + locs.resize(points.size()); + for (size_t i = 0; i < points.size(); i++) + if ((locs[i].dist = this->findPoint(points[i],locs[i].u)) >= 0.0) + ++nFound; + + return nFound == points.size(); +} + + bool ASMbase::getSolution (Matrix& sField, const Vector& locSol, const IntVec& nodes) const { @@ -1772,7 +1785,7 @@ bool ASMbase::evalSolution (Matrix&, const IntegrandBase&, bool ASMbase::evalSolution (Matrix& sField, const IntegrandBase& integrand, - const IntVec& elements) const + const IntVec& elements, const RealArray*) const { if (elements.empty()) return this->evalSolution(sField,integrand,nullptr,true); diff --git a/src/ASM/ASMbase.h b/src/ASM/ASMbase.h index 6aec0d545..ab3687b97 100644 --- a/src/ASM/ASMbase.h +++ b/src/ASM/ASMbase.h @@ -317,6 +317,8 @@ class ASMbase void printNodes(std::ostream& os) const; //! \brief Prints out element connections of this patch to the given stream. void printElements(std::ostream& os) const; + //! \brief Prints out additional app-dependent element information. + virtual void printElmInfo(int, const IntegrandBase*) const {} //! \brief Increase all global node numbers by \a nshift. virtual void shiftGlobalNodeNums(int nshift); @@ -582,6 +584,23 @@ class ASMbase //! \return Distance from the point \a X to the found point virtual double findPoint(Vec3& X, double* param) const = 0; + //! \brief A struct with element and associated local coordinates of a point. + struct PointParams + { + size_t iel; //!< 1-based element index + double u[3]; //!< Parametric coordinates w.r.t. element or patch + double dist; //!< Distance to spatial point + + //! \brief Default constructor. + explicit PointParams(size_t e = 0) : iel(e), u{0.0,0.0,0.0}, dist(-1.0) {} + }; + + //! \brief Finds elements and local coordinates matching spatial points. + //! \param points List of spatial point coordinates + //! \param[out] locs List of matching elements and/or parametric coordinates + virtual bool findPoints(std::vector& points, + std::vector& locs) const; + //! \brief Creates a standard FE model of this patch for visualization. //! \param[out] grid The generated finite element grid //! \param[in] npe Number of visualization nodes over each knot span @@ -714,13 +733,17 @@ class ASMbase //! \param[out] sField Solution field //! \param[in] integrand Object with problem-specific data and methods //! \param[in] elements List of elements to evaluate for (all if empty) + //! \param[in] lpar Local parameter values of the result sampling points //! //! \details The secondary solution is derived from the primary solution, //! which is assumed to be stored within the \a integrand for current patch. //! This method evaluates the secondary solution at the center of the //! provided list of \a elements, or for all if the list is empty. + //! If \a lpar is specified, it is assumed to contain the local parameters + //! (w.r.t. the elements) of the evaluation points virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand, - const IntVec& elements) const; + const IntVec& elements, + const RealArray* lpar = nullptr) const; //! \brief Projects the secondary solution using a (discrete) global L2-fit. //! \param[out] sField Secondary solution field control point values diff --git a/src/ASM/ASMs2DLag.C b/src/ASM/ASMs2DLag.C index 0fb25a6d4..4d7a3f045 100644 --- a/src/ASM/ASMs2DLag.C +++ b/src/ASM/ASMs2DLag.C @@ -718,7 +718,8 @@ bool ASMs2DLag::tesselate (ElementBlock& grid, const int*) const } -void ASMs2DLag::constrainEdge (int dir, bool open, int dof, int code, char basis) +void ASMs2DLag::constrainEdge (int dir, bool open, int dof, int code, + char basis) { this->ASMs2D::constrainEdge(dir, open, dof, code > 0 ? -code : code, basis); } @@ -741,39 +742,79 @@ bool ASMs2DLag::evalSolution (Matrix& sField, const Vector& locSol, size_t nCmp = locSol.size() / this->getNoProjectionNodes(); size_t ni = gpar ? gpar[0].size() : nel; size_t nj = gpar ? gpar[1].size() : 1; - size_t nen = p1*p2; sField.resize(nCmp,ni*nj); - Matrix elmSol(nCmp,nen); - RealArray N(nen), val; + RealArray N, val; double xi = 0.0, eta = 0.0; if (!gpar && !Lagrange::computeBasis(N,p1,xi,p2,eta)) return false; - size_t ip = 1; + size_t ip = 0; int iel = 0; for (size_t j = 0; j < nj; j++) - for (size_t i = 0; i < ni; i++, ip++) + for (size_t i = 0; i < ni; i++) { if (gpar) - { iel = this->findElement(gpar[0][i], gpar[1][j], &xi, &eta); - if (iel < 1 || iel > static_cast(nel)) - return false; - if (!Lagrange::computeBasis(N,p1,xi,p2,eta)) - return false; - } - else - iel++; + else // evaluate at element centers + --iel; - for (size_t a = 1; a <= nen; a++) - elmSol.fillColumn(a, locSol.ptr() + nCmp*MNPC[iel-1][a-1]); - - elmSol.multiply(N,val); - sField.fillColumn(ip,val); + if (!this->evalSolPt(iel,xi,eta,nCmp,locSol,val,N)) + return false; + else if (!val.empty()) // skip elements with no results + sField.fillColumn(++ip,val); } + sField.resize(nCmp,ip); + return true; +} + + +bool ASMs2DLag::evalSolution (Matrix& sField, const Vector& locSol, + const IntVec& elms, const RealArray* lpar) const +{ + const size_t nCmp = locSol.size() / this->getNoProjectionNodes(); + const size_t nPts = lpar ? lpar[0].size() : elms.size(); + const size_t iPts = sField.cols(); + sField.resize(nCmp,iPts+nPts); + + RealArray N, ptSol; + for (size_t i = 0; i < nPts; i++) + { + double xi = lpar ? lpar[0][i] : 0.0; + double eta = lpar ? lpar[1][i] : 0.0; + if (!this->evalSolPt(elms[i],xi,eta,nCmp,locSol,ptSol,N)) + return false; + + sField.fillColumn(iPts+i+1,ptSol); + } + + return true; +} + + +bool ASMs2DLag::evalSolPt (int iel, double xi, double eta, size_t nCmp, + const Vector& pchSol, RealArray& ptSol, + RealArray& N) const +{ + if (iel < 0 && -iel <= static_cast(nel)) // assume N is already evaluated + iel = -iel-1; + else if (iel < 1 || iel > static_cast(nel)) + return false; + else if (!Lagrange::computeBasis(N,p1,xi,p2,eta)) + return false; + else + --iel; + + ptSol.assign(nCmp,0.0); + for (size_t a = 0; a < N.size(); a++) + { + size_t ip = nCmp*MNPC[iel][a]; + for (size_t i = 0; i < nCmp; i++) + ptSol[i] += N[a]*pchSol[ip++]; + } + return true; } @@ -997,7 +1038,7 @@ int ASMs2DLag::findElement (double u, double v, double* xi, double* eta) const if (!surf) { std::cerr <<" *** ASMs2DLag::findElement: No spline geometry"<< std::endl; - return -1; + return 0; } int ku, kv; diff --git a/src/ASM/ASMs2DLag.h b/src/ASM/ASMs2DLag.h index a25229e47..76b9c7d70 100644 --- a/src/ASM/ASMs2DLag.h +++ b/src/ASM/ASMs2DLag.h @@ -233,6 +233,20 @@ class ASMs2DLag : public ASMs2D, protected ASMLagBase const RealArray* gpar, bool regular, int, int) const; + //! \brief Evaluates the primary solution field at the given points. + //! \param sField Solution field + //! \param[in] locSol Solution vector local to current patch + //! \param[in] elms List of elements to evaluate at + //! \param[in] lpar Local parameter values of the result sampling points + //! + //! \details If \a lpar is null, the solution is evaluated at the center of + //! each element. Otherwise, it is assumed to contain the local ξ and η + //! parameters of the evaluation points, w.r.t. the provided elements. + //! \note The resuls are appended as extra columns to output matrix \a sField. + //! Therefore, any existing results is preserved. + bool evalSolution(Matrix& sField, const Vector& locSol, + const IntVec& elms, const RealArray* lpar = nullptr) const; + //! \brief Evaluates the secondary solution field at all visualization points. //! \param[out] sField Solution field //! \param[in] integrand Object with problem-specific data and methods @@ -290,6 +304,11 @@ class ASMs2DLag : public ASMs2D, protected ASMLagBase static void createMNPC(size_t nx, size_t ny, int p1, int p2, IntMat& MNPC); protected: + //! \brief Evaluates a nodal solution field at specified point in an element. + virtual bool evalSolPt(int iel, double xi, double eta, size_t nCmp, + const Vector& pchSol, RealArray& ptSol, + RealArray& N) const; + size_t nx; //!< Number of nodes in first parameter direction size_t ny; //!< Number of nodes in second parameter direction int p1; //!< Polynomial order in first parameter direction diff --git a/src/ASM/ASMs3DLag.C b/src/ASM/ASMs3DLag.C index c210bf8c0..4809e0c80 100644 --- a/src/ASM/ASMs3DLag.C +++ b/src/ASM/ASMs3DLag.C @@ -952,7 +952,7 @@ int ASMs3DLag::findElement (double u, double v, double w, if (!svol) { std::cerr <<" *** ASMs3DLag::findElement: No spline geometry"<< std::endl; - return -1; + return 0; } std::array,3> knot; @@ -1001,11 +1001,9 @@ bool ASMs3DLag::evalSolution (Matrix& sField, const Vector& locSol, size_t ni = gpar ? gpar[0].size() : nel; size_t nj = gpar ? gpar[1].size() : 1; size_t nk = gpar ? gpar[2].size() : 1; - size_t nen = p1*p2*p3; sField.resize(nCmp,ni*nj*nk); - Matrix elmSol(nCmp,nen); - RealArray N(nen), val; + RealArray N, val; double xi = 0.0, eta = 0.0, zeta = 0.0; if (!gpar && !Lagrange::computeBasis(N,p1,xi,p2,eta,p3,zeta)) @@ -1018,21 +1016,13 @@ bool ASMs3DLag::evalSolution (Matrix& sField, const Vector& locSol, for (size_t i = 0; i < ni; i++, ip++) { if (gpar) - { iel = this->findElement(gpar[0][i], gpar[1][j], gpar[2][k], &xi, &eta, &zeta); - if (iel < 1 || iel > static_cast(nel)) - return false; - if (!Lagrange::computeBasis(N,p1,xi,p2,eta,p3,zeta)) - return false; - } - else - iel++; - - for (size_t a = 1; a <= nen; a++) - elmSol.fillColumn(a, locSol.ptr() + nCmp*MNPC[iel-1][a-1]); + else // evaluate at element centers + --iel; + if (!this->evalSolPt(iel,xi,eta,zeta,nCmp,locSol,val,N)) + return false; - elmSol.multiply(N,val); sField.fillColumn(ip,val); } @@ -1040,6 +1030,31 @@ bool ASMs3DLag::evalSolution (Matrix& sField, const Vector& locSol, } +bool ASMs3DLag::evalSolPt (int iel, double xi, double eta, double zeta, + size_t nCmp, const Vector& pchSol, RealArray& ptSol, + RealArray& N) const +{ + if (iel < 0 && -iel <= static_cast(nel)) // assume N is already evaluated + iel = -iel-1; + else if (iel < 1 || iel > static_cast(nel)) + return false; + else if (!Lagrange::computeBasis(N,p1,xi,p2,eta,p3,zeta)) + return false; + else + --iel; + + ptSol.assign(nCmp,0.0); + for (size_t a = 0; a < N.size(); a++) + { + size_t ip = nCmp*MNPC[iel][a]; + for (size_t i = 0; i < nCmp; i++) + ptSol[i] += N[a]*pchSol[ip++]; + } + + return true; +} + + bool ASMs3DLag::evalSolution (Matrix& sField, const IntegrandBase& integrand, const int*, char) const { diff --git a/src/ASM/ASMs3DLag.h b/src/ASM/ASMs3DLag.h index c49e5695e..f0ae72348 100644 --- a/src/ASM/ASMs3DLag.h +++ b/src/ASM/ASMs3DLag.h @@ -312,6 +312,10 @@ class ASMs3DLag : public ASMs3D, protected ASMLagBase static void createMNPC(size_t nx, size_t ny, size_t nz, int p1, int p2, int p3, IntMat& MNPC); + //! \brief Evaluates a nodal solution field at specified point in an element. + bool evalSolPt(int iel, double xi, double eta, double zeta, size_t nCmp, + const Vector& pchSol, RealArray& ptSol, RealArray& N) const; + protected: size_t nx; //!< Number of nodes in first parameter direction size_t ny; //!< Number of nodes in second parameter direction diff --git a/src/ASM/ASMu2DLag.C b/src/ASM/ASMu2DLag.C index 4d5509561..eb5e079ad 100644 --- a/src/ASM/ASMu2DLag.C +++ b/src/ASM/ASMu2DLag.C @@ -21,11 +21,30 @@ #include "IFEM.h" #include #include +#include #ifdef USE_OPENMP #include #endif +extern "C" { + //! \brief Parametric inversion for a linear triangle. + void t3_inv_(const int& nsd, const double& tol, const double& eps, + const double* X, const double* Xn, double* xi, + const int& ipsw, const int& iwr, int& ierr); + //! \brief Parametric inversion for a bilinear quadrilateral. + void q4_inv_(const int& nsd, const double& tol, const double& eps, + const double* X, const double* Xn, double* xi, + const int& ipsw, const int& iwr, int& ierr); + //! \brief Evaluate the bilinear shape functions. + void shape2_(const int& iop, const int& nenod, const double* xi, + double* shp, const int& ipsw, const int& iwr, int& ierr); + //! \brief Opens a temporary file for Fortran print. + int openftnfile_(const char* fname, const ssize_t nchar); + //! \brief Closes a Fortran file. + void closeftnfile_(const int& iunit); +} + ASMu2DLag::ASMu2DLag (unsigned char n_s, unsigned char n_f, char fType) : ASMs2DLag(n_s,n_f) @@ -587,3 +606,229 @@ bool ASMu2DLag::integrate (Integrand& integrand, return ok; } + + +bool ASMu2DLag::evalSolPt (int iel, double xi, double eta, size_t nCmp, + const Vector& pchSol, RealArray& ptSol, + RealArray& N) const +{ + ptSol.clear(); + size_t jel = iel > 0 ? iel : -iel; + if (jel < 1 || jel > MNPC.size()) + return false; + else if (MNPC[--jel].size() > 3) // 4-noded 2D element + return this->ASMs2DLag::evalSolPt(iel,xi,eta,nCmp,pchSol,ptSol,N); + else if (MNPC[jel].size() < 3) + return true; // no results for 0D and 1D elements + + // This element is a linear triangle + RealArray N3(3,1.0/3.0); + if (iel > 0) N3 = { xi, eta, 1.0-xi-eta }; + return this->ASMs2DLag::evalSolPt(-1-jel,xi,eta,nCmp,pchSol,ptSol,N3); +} + + +bool ASMu2DLag::writeXML (const char* fname) const +{ + std::ofstream os(fname); + if (!os) return false; + + os <<"\n "; + for (const Vec3& X : coord) + os <<"\n "<< X; + os <<"\n "; + for (size_t nen = 1; nen <= 4; nen++) + { + bool haveElms = false; + for (const IntVec& mnpc : MNPC) + if (mnpc.size() == nen) + { + if (!haveElms) + os <<"\n "; + os <<"\n "; + for (int inod : mnpc) + os <<" "<< inod; + haveElms = true; + } + if (haveElms) + os <<"\n "; + } + os <<"\n\n"; + + return true; +} + + +bool ASMu2DLag::findPoints (Vec3Vec& points, ParamsVec& locs) const +{ + // Calculate the bounding box for all (shell) elements + size_t nonshell = 0; + std::vector bbox(nel); + for (size_t iel = 0; iel < nel; iel++) + if (MNPC[iel].size() == 3 || MNPC[iel].size() == 4) + this->getBoundingBox(MNPC[iel],bbox[iel]); + else + { + nonshell++; // make sure any other elements are sorted last + bbox[iel].first = bbox[iel].second = Vec3(1.0e99,1.0e99,1.0e99); + } + + std::array,3> eIdx; + std::vector::iterator it; + + // Index sort the elements in each spatial direction w.r.t. the max coordinate + for (int i = 0; i < 3; i++) + { + eIdx[i].resize(nel); + std::iota(eIdx[i].begin(), eIdx[i].end(), 0); + std::sort(eIdx[i].begin(), eIdx[i].end(), [i,&bbox](size_t a, size_t b) + { return bbox[a].second[i] < bbox[b].second[i]; }); + if (nonshell > 0) + eIdx[i].resize(nel-nonshell); // Remove the non-shell elements +#if SP_DEBUG > 2 + std::cout <<"\nElement order in "<< char('X'+i) <<"-direction:"; + for (size_t iel = 0; iel < nel-nonshell; iel++) + std::cout << (iel%10 ? ' ' : '\n') << MLGE[eIdx[i][iel]] + <<" ("<< bbox[eIdx[i][iel]].second[i] <<")"; + std::cout << std::endl; +#endif + } + + Vec3Pair BBox; // Bounding box for the whole patch + double maxDim = 0.0; + for (int i = 0; i < 3; i++) + { + BBox.first[i] = bbox[eIdx[i].front()].first[i]; + BBox.second[i] = bbox[eIdx[i].back()].second[i]; + if (double dim = BBox.second[i] - BBox.first[i]; dim > maxDim) + maxDim = dim; + } +#if SP_DEBUG > 1 + std::cout <<"\nDomain bounding box: [ " + << BBox.first <<" , "<< BBox.second + <<" ] maxDim = "<< maxDim << std::endl; +#endif + const double tol = maxDim*1.0e-8; + BBox.second += Vec3(RealArray(3,-tol)); + + size_t nFound = 0; + locs.clear(); + locs.reserve(points.size()); + for (Vec3& X : points) + { + // Do a binary search in each coordinate direction + // to find the range of elements to check in more detail + std::set nearElms; + for (int i = 0; i < 3; i++) + if (BBox.second[i] > BBox.first[i]) // skip the flat dimension + { + it = std::upper_bound(eIdx[i].begin(), eIdx[i].end(), X[i], + [i,&bbox](double x, size_t a) + { return x < bbox[a].second[i]; }); + if (nearElms.empty()) + nearElms.insert(it,eIdx[i].end()); + else + { + std::set tmp1(it,eIdx[i].end()), tmp2; + nearElms.swap(tmp2); + std::set_intersection(tmp1.begin(), tmp1.end(), + tmp2.begin(), tmp2.end(), + std::inserter(nearElms,nearElms.begin())); + } + } + + // Lambda function doing the Y < X comparison for two Vec3 objects + auto&& below = [&X,tol](const Vec3& Y) + { + for (int i = 0; i < 3; i++) + if (Y[i] > X[i]+tol) + return false; + return true; + }; + +#if SP_DEBUG > 1 + std::cout <<"\nSearching spatial point "<< X <<" (max "<< nearElms.size() + <<" elements to check)."<< std::endl; +#endif + bool found = false; + for (size_t iel : nearElms) + if (below(bbox[iel].first)) + { +#if SP_DEBUG > 1 + std::cout <<"Checking element "<< MLGE[iel] <<" in [ "<< bbox[iel].first + <<" - "<< bbox[iel].second <<" ]"<< std::endl; +#endif + // Perform parametric inversion to see if this point + // really is within current element + if (PointParams elem(1+iel); this->paramInvert(X,elem)) + { + locs.emplace_back(elem); + found = true; + break; + } + } + + if (found) + nFound++; + else + { + locs.emplace_back(PointParams()); + IFEM::cout <<" ** Warning: No element matches the point "<< X + << std::endl; + } + } + + return nFound == points.size(); +} + + +bool ASMu2DLag::paramInvert (Vec3& X, PointParams& elm) const +{ + Matrix Xn; + if (!this->getElementCoordinates(Xn,elm.iel)) + return false; + + const int nelnod = Xn.cols(); + if (nelnod == 4) // swap element node 3 and 4 + for (size_t r = 1; r <= Xn.rows(); r++) + std::swap(Xn(r,3),Xn(r,4)); + + const double eps = 1.0e-15; + const double tol = 1.0e-3; +#if SP_DEBUG > 2 + int ipsw = SP_DEBUG; +#else + int ipsw = 0; +#endif + // Open a temporary file for the Fortran output + std::string tmpfile = "/tmp/" + std::string(getenv("USER")) + ".ftn"; + const int iwr = openftnfile_(tmpfile.c_str(),tmpfile.size()); + if (iwr <= 0) ipsw = -1; // suppress all output + int ierr = 0; + if (nelnod == 3) + t3_inv_(nsd,tol,eps,X.ptr(),Xn.ptr(),elm.u,ipsw,iwr,ierr); + else if (nelnod == 4) + q4_inv_(nsd,tol,eps,X.ptr(),Xn.ptr(),elm.u,ipsw,iwr,ierr); + else + ierr = -99; + if (iwr > 0 && ipsw >= 0) + { + // Copy fortran output to IFEM::cout + closeftnfile_(iwr); + std::string cline; + std::ifstream is(tmpfile); + while (std::getline(is,cline)) + IFEM::cout << cline << std::endl; + } + if (ierr != 0 && ierr != 2) + return false; + + // Calculate the distance to the projected point + Vector N(nelnod); + shape2_(0,nelnod,elm.u,N.ptr(),ipsw,iwr,ierr); + Vec3 X0(Xn * N); + elm.dist = (X0-X).length(); + + X = X0; // update the spatial coordinates + return true; +} diff --git a/src/ASM/ASMu2DLag.h b/src/ASM/ASMu2DLag.h index 1e2f5be23..69cca9f4e 100644 --- a/src/ASM/ASMu2DLag.h +++ b/src/ASM/ASMu2DLag.h @@ -136,10 +136,28 @@ class ASMu2DLag : public ASMs2DLag //! \note The number of element nodes must be set in \a grid on input. virtual bool tesselate(ElementBlock& grid, const int*) const; + using ParamsVec = std::vector; //!< Convenience type alias + + //! \brief Finds elements and local coordinates matching spatial points. + //! \param points List of spatial point coordinates + //! \param[out] locs List of matching element and local coordinates + virtual bool findPoints(Vec3Vec& points, ParamsVec& locs) const; + + //! \brief Dumps the mesh to the specified XML-file. + bool writeXML(const char* fname) const; + protected: + //! \brief Evaluates a nodal solution field at specified point in an element. + virtual bool evalSolPt(int iel, double xi, double eta, size_t nCmp, + const Vector& pchSol, RealArray& ptSol, + RealArray& N) const; + //! \brief Generate thread groups using multi-coloring. void generateThreadGroupsMultiColored(bool silence, bool separateGroup1Noded); + //! \brief Performs parametric inversion for an element. + bool paramInvert(Vec3& X, PointParams& elm) const; + bool swapNode34; //!< If \e true, element nodes 3 and 4 should be swapped std::vector nodeSets; //!< Node sets for Dirichlet BCs diff --git a/src/ASM/Ftn/dcros3.f b/src/ASM/Ftn/dcros3.f new file mode 100644 index 000000000..2cc69e287 --- /dev/null +++ b/src/ASM/Ftn/dcros3.f @@ -0,0 +1,12 @@ + SUBROUTINE DCROS3 (X,Y,Z) +C +C Calculate {Z} = {X} x {Y} +C + DOUBLE PRECISION X(3),Y(3),Z(3) +C + Z(1) = X(2)*Y(3) - Y(2)*X(3) + Z(2) = X(3)*Y(1) - Y(3)*X(1) + Z(3) = X(1)*Y(2) - Y(1)*X(2) +C + RETURN + END diff --git a/src/ASM/Ftn/dmuax3.f b/src/ASM/Ftn/dmuax3.f new file mode 100644 index 000000000..582f6f1fc --- /dev/null +++ b/src/ASM/Ftn/dmuax3.f @@ -0,0 +1,12 @@ + SUBROUTINE DMUAX3 (A,X,Y) +C +C Calculate {Y} = [A] * {X} +C + DOUBLE PRECISION A(3,3),X(3),Y(3) +C + Y(1) = A(1,1)*X(1) + A(1,2)*X(2) + A(1,3)*X(3) + Y(2) = A(2,1)*X(1) + A(2,2)*X(2) + A(2,3)*X(3) + Y(3) = A(3,1)*X(1) + A(3,2)*X(2) + A(3,3)*X(3) +C + RETURN + END diff --git a/src/ASM/Ftn/dnorm3.f b/src/ASM/Ftn/dnorm3.f new file mode 100644 index 000000000..55c5fa2da --- /dev/null +++ b/src/ASM/Ftn/dnorm3.f @@ -0,0 +1,33 @@ + DOUBLE PRECISION FUNCTION DNORM3 (V) + DOUBLE PRECISION V(3),DV +C +C DNORM3 : Calculate |v| and 1/|v| * {v} +C + DOUBLE PRECISION ZERO , EPS , ONE + PARAMETER ( ZERO = 0.0D0, EPS = 1.0D-15, ONE = 1.0D0 ) +C + DV = SQRT(V(1)*V(1) + V(2)*V(2) + V(3)*V(3)) + IF (DV .LT. EPS) THEN + DNORM3 = ZERO + ELSE + DNORM3 = DV + V(1) = V(1) / DV + V(2) = V(2) / DV + V(3) = V(3) / DV + IF (ABS(V(1)) .GT. ONE-EPS) THEN + V(1) = SIGN(ONE,V(1)) + V(2) = ZERO + V(3) = ZERO + ELSE IF (ABS(V(2)) .GT. ONE-EPS) THEN + V(1) = ZERO + V(2) = SIGN(ONE,V(2)) + V(3) = ZERO + ELSE IF (ABS(V(3)) .GT. ONE-EPS) THEN + V(1) = ZERO + V(2) = ZERO + V(3) = SIGN(ONE,V(3)) + ENDIF + ENDIF +C + RETURN + END diff --git a/src/ASM/Ftn/dprint.f b/src/ASM/Ftn/dprint.f new file mode 100644 index 000000000..383baae62 --- /dev/null +++ b/src/ASM/Ftn/dprint.f @@ -0,0 +1,49 @@ + SUBROUTINE DPRINT (A,M,N,CHEAD,IWR) + INTEGER M,N, IWR + DOUBLE PRECISION A(M,*) + CHARACTER CHEAD*(*) +C +C PURPOSE: +C DPRINT : Print the double precision array A(M,N). +C +C METHOD: +C The array is printed as a two-dimensional matrix. +C Rows and columns are numbered and each element is printed with +C the format 1PE14.6, i.e. 7 significant digits. +C A heading (e.g. the array name) is also printed. +C +C INPUT ARGUMENTS: +C A - Array to be printed +C M - Number of rows +C N - Number of columns +C CHEAD - Heading +C IWR - Print unit number +C + INTEGER NCOL,J,I,K,L + PARAMETER ( NCOL = 9 ) +C + IF (N .LE. 0 .OR. M .LE. 0) GOTO 8000 +C + WRITE(IWR,9000) CHEAD +C + DO 2000 J = 1,N,NCOL + K = MIN(N,J+NCOL-1) +C + WRITE(IWR,9100) (L,L=J,K) + WRITE(IWR,9200) ('--------------',L=J,K) +C + DO 1000 I = 1,M + WRITE(IWR,9300) I,(A(I,L),L=J,K) + 1000 CONTINUE +C + 2000 CONTINUE +C + 8000 CONTINUE + RETURN +C + 9000 FORMAT( / 5X, A, ' = ' ) + 9100 FORMAT( / 4X, 9(I9,5X) ) + 9200 FORMAT( 4X, 9A14 ) + 9300 FORMAT( I3, '|', 1P,9E14.6 ) +C + END diff --git a/src/ASM/dslbln.f b/src/ASM/Ftn/dslbln.f similarity index 100% rename from src/ASM/dslbln.f rename to src/ASM/Ftn/dslbln.f diff --git a/src/ASM/dsolv3.f b/src/ASM/Ftn/dsolv3.f similarity index 100% rename from src/ASM/dsolv3.f rename to src/ASM/Ftn/dsolv3.f diff --git a/src/ASM/Ftn/errmsg.f b/src/ASM/Ftn/errmsg.f new file mode 100644 index 000000000..9851b0dea --- /dev/null +++ b/src/ASM/Ftn/errmsg.f @@ -0,0 +1,42 @@ + SUBROUTINE ERRMSG (PRGNAM,IWR,IERR) + CHARACTER PRGNAM*(*) + INTEGER IWR,IERR +C +C ERRMSG : Print an error/warning message. +C +C PRGNAM - Name of subroutine where the error occurred +C IERR - Error flag +C > 1 : Warning +C = -1 : Traceback +C < -1 : Error +C + IF (IERR .LT. -1) THEN + WRITE(IWR,9901) PRGNAM + ELSE IF (IERR .GT. 1) THEN + WRITE(IWR,9902) PRGNAM + ELSE IF (IERR .EQ. -1) THEN + WRITE(IWR,9903) PRGNAM + ENDIF +C + IF (ABS(IERR) .LE. 1) GOTO 8000 +C + IF (IERR .EQ. -2) THEN + WRITE(IWR,9904) 'INVALID INPUT ARGUMENTS' + ELSE IF (IERR .EQ. -98) THEN + WRITE(IWR,9904) 'ARRAY IS NOT ALLOCATED' + ELSE IF (IERR .EQ. -99) THEN + WRITE(IWR,9904) 'UNABLE TO ALLOCATE MEMORY' + ELSE + WRITE(IWR,9905) IERR + ENDIF +C + 8000 CONTINUE + RETURN +C + 9901 FORMAT( / ' *** ERROR RETURN FROM SUBROUTINE ', A ) + 9902 FORMAT( / ' ** WARNING FROM SUBROUTINE ', A ) + 9903 FORMAT( ' * TRACEBACK: ', A ) + 9904 FORMAT( 5X, A ) + 9905 FORMAT( 5X, 'ERROR FLAG = ', I6 ) +C + END diff --git a/src/ASM/Ftn/q4_inv.f b/src/ASM/Ftn/q4_inv.f new file mode 100644 index 000000000..ee46ebb02 --- /dev/null +++ b/src/ASM/Ftn/q4_inv.f @@ -0,0 +1,274 @@ + SUBROUTINE Q4_INV (NDIM,EPSO,EPSZ,X,TENC,XI,IPSW,IWR,IERR) +C +C $Id$ +C***************************** USERS' PART ***************************** +C +C A F E M ROUTINE : Q4_INV +C +C PURPOSE: +C Q4_INV : Find natural coordinates for a point within a Q4-element. +C +C METHOD: +C The natural coordinates {XI} corresponding to the +C point {X} are calculated. The element nodes are assumed to have +C natural coordinates (-1,-1), (1,-1), (1,1), (-1,1), respectively. +C +C INPUT ARGUMENTS: +C NDIM - Number of spatial dimensions +C EPSO - Tolerence for checking if point is outside the element +C EPSZ - Tolerance for equality check with zero +C X - Cartesian coordinates of a point within the element +C TENC - Table of element nodal coordinates +C IPSW - Print switch +C IWR - Print unit number +C +C INPUT/OUTPUT ARGUMENTS: +C None +C +C OUTPUT ARGUMENTS: +C XI - Natural coordinates of the point {X} +C IERR - Error flag +C +C RESTRICTIONS: +C None +C +C************************** MAINTENANCE PART *************************** +C +C INTERNAL VARIABLES: +C None +C +C PRINT SWITCH: +C .LT. 0 - No print +C .EQ. 0 - Print error messages only +C .GE. 7 - Print ENTER and LEAVE +C .GE. 8 - Print in addition INPUT parameters +C .GE. 9 - Print in addition OUTPUT parameters +C +C SUBROUTINES CALLED: +C DDOT (double precision function) +C DMUAX3 +C DSLBLN +C SHAPE2 +C SHLTRI +C +C CODING: +C 19-SEP-96, Knut Morten Okstad, NTNU +C +C REVISION: +C +C********************** END SUBROUTINE DESCRIPTION ********************* +C + IMPLICIT NONE +C +C * Global variables +C + INTEGER NDIM,IPSW,IWR,IERR + DOUBLE PRECISION EPSO,EPSZ,X(NDIM),TENC(NDIM,4),XI(2) +C +C * Internal variables + LOGICAL LXY + INTEGER i,IERRL,j,NSOL + DOUBLE PRECISION A(4,2),D1,D2,TOL,TGL(9),TLX(2,5), + + X1(4),X2(4),XL(3) +C + DOUBLE PRECISION ZERO , C , ONE + PARAMETER ( ZERO = 0.0D0, C = 0.25D0, ONE = 1.0D0 ) +C + DOUBLE PRECISION DDOT +C +C**************************** END VARIABLES **************************** +C +C Entry section +C + IF (IPSW .GE. 7) WRITE(IWR,9900) + IF (IPSW .GE. 8) THEN + WRITE(IWR,9920) EPSO,EPSZ,X + CALL DPRINT (TENC,NDIM,4,'TENC ',IWR) + ENDIF +C +C +C Program logic section +C +C * Check input arguments + IERR = -1 + IF (EPSZ .LE. ZERO .OR. EPSO .LT. EPSZ) GOTO 7200 +C + IF (NDIM .EQ. 2) THEN + LXY = .TRUE. + ELSE IF (NDIM .EQ. 3) THEN + TOL = MAX(TENC(3,1),TENC(3,2),TENC(3,3),TENC(3,4)) + + - MIN(TENC(3,1),TENC(3,2),TENC(3,3),TENC(3,4)) + LXY = TOL .LT. EPSZ + ELSE + GOTO 7200 + ENDIF +C + IF (LXY) THEN +C +C * The element is parallel +C to the global XY-plane + DO 1000 j = 1,2 + A(1,j) = ( TENC(j,1) - TENC(j,2) + TENC(j,3) - TENC(j,4))*C + A(2,j) = (-TENC(j,1) + TENC(j,2) + TENC(j,3) - TENC(j,4))*C + A(3,j) = (-TENC(j,1) - TENC(j,2) + TENC(j,3) + TENC(j,4))*C + A(4,j) = ( TENC(j,1) + TENC(j,2) + TENC(j,3) + TENC(j,4))*C + 1000 CONTINUE +C + A(4,1) = X(1) - A(4,1) + A(4,2) = X(2) - A(4,2) +C + ELSE +C +C * The element is not parallel +C to the global XY-plane, need +C to calculate local coordinates + DO 2000 j = 1,NDIM + X1(j) = TENC(j,2) - TENC(j,1) + X2(j) = TENC(j,3) - TENC(j,1) + 2000 CONTINUE + CALL SHLTRI (2,X1,X2,TGL,IPSW,IWR,IERRL) + IF (IERRL .LT. 0) GOTO 7100 +C + DO 2100 j = 1,4 + CALL DMUAX3 (TGL,TENC(1,j),TLX(1,j)) + 2100 CONTINUE + CALL DMUAX3 (TGL,X,XL) +C +C * Initiate the aux. variables + DO 2200 j = 1,2 + A(1,j) = ( TLX(j,1) - TLX(j,2) + TLX(j,3) - TLX(j,4)) * C + A(2,j) = (-TLX(j,1) + TLX(j,2) + TLX(j,3) - TLX(j,4)) * C + A(3,j) = (-TLX(j,1) - TLX(j,2) + TLX(j,3) + TLX(j,4)) * C + A(4,j) = ( TLX(j,1) + TLX(j,2) + TLX(j,3) + TLX(j,4)) * C + 2200 CONTINUE +C + A(4,1) = XL(1) - A(4,1) + A(4,2) = XL(2) - A(4,2) +C + ENDIF +C +C----------------------------------------------------------------------- +C +C We have to solve the following set of equations (j=1,2): +C +C A1j*XI*ETA + A2j*XI + A3j*ETA = A4j +C +C The way that we may solve this nonlinear (in XI and ETA) +C set of equations depends on the coefficients Aij. +C The solution is unique for proper input. See DSLBLN for details. +C +C----------------------------------------------------------------------- +C + CALL DSLBLN (IPSW,IWR,EPSZ,A(1,1),A(1,2),NSOL,X1,X2) + IF (NSOL .LE. 0) GOTO 7100 +C +C * Check that the solutions +C are "inside" the element + TOL = ONE + EPSO + DO 3000 j = 1,NSOL + IF (ABS(X1(j)).LT.TOL .AND. ABS(X2(j)).LT.TOL) THEN + IERR = IERR + 1 + IF (IERR .GT. 0) THEN +C +C * We have multiple solutions. +C Check if they are almost equal. +C + IF (ABS(XI(1)-X1(j)) .LT. EPSZ*(TOL+TOL) .AND. + + ABS(XI(2)-X2(j)) .LT. EPSZ*(TOL+TOL)) THEN + XI(1) = 0.5D0*(XI(1)+X1(j)) + XI(2) = 0.5D0*(XI(2)+X2(j)) + IERR = IERR - 1 + GOTO 3000 + ENDIF +C +C * Choose the one closest to X +C + CALL SHAPE2 (0,4,XI,TGL,IPSW,IWR,IERRL) + IF (IERRL .LT. 0) GOTO 7100 +C + D1 = ZERO + DO 2500 i = 1,NDIM + D1 = D1 + (X(i)-DDOT(4,TGL,1,TENC(i,1),NDIM))**2.0D0 + 2500 CONTINUE +C + XL(1) = X1(j) + XL(2) = X2(j) + CALL SHAPE2 (0,4,XL,TGL,IPSW,IWR,IERRL) + IF (IERRL .LT. 0) GOTO 7100 +C + D2 = ZERO + DO 2600 i = 1,NDIM + D2 = D2 + (X(i)-DDOT(4,TGL,1,TENC(i,1),NDIM))**2.0D0 + 2600 CONTINUE + IF (D1 .LT. D2) GOTO 3000 +C + ENDIF + XI(1) = X1(j) + XI(2) = X2(j) + ENDIF + 3000 CONTINUE + + IF (IERR .GT. 0) GOTO 6200 + IF (IERR .EQ. 0) GOTO 8000 +C +C +C Warning section +C + 6100 CONTINUE + IERR = 1 + IF (IPSW .GE. 0) THEN + WRITE(IWR,9961) (j,X1(j),X2(j),j=1,NSOL) + WRITE(IWR,9920) EPSO,EPSZ,X + CALL DPRINT (TENC,NDIM,4,'TENC ',IWR) + ENDIF + GOTO 8000 + 6200 CONTINUE + IERR = 2 + IF (IPSW .GE. 0) THEN + WRITE(IWR,9962) (j,X1(j),X2(j),j=1,NSOL) + WRITE(IWR,9963) XI,MIN(D1,D2) + WRITE(IWR,9920) EPSO,EPSZ,X + CALL DPRINT (TENC,NDIM,4,'TENC ',IWR) + ENDIF + GOTO 8000 +C +C Error section +C + 7100 CONTINUE + IERR = -1 + GOTO 7900 + 7200 CONTINUE + IERR = -2 +C + 7900 CONTINUE + IF (IPSW .GE. 0) THEN + CALL ERRMSG ('Q4_INV',IWR,IERR) + ENDIF +C +C +C Closing section +C + 8000 CONTINUE + IF (IPSW .GE. 7) WRITE(IWR,9910) + IF (IPSW .GE. 9) WRITE(IWR,9930) IERR,XI +C + RETURN +C + 9900 FORMAT( / ' ENTERING SUBROUTINE Q4_INV' ) + 9910 FORMAT( / ' LEAVING SUBROUTINE Q4_INV' ) + 9920 FORMAT( ' WITH INPUT PARAMETERS:' + + / ' EPS = ', 1P,2E13.5, + + / ' X = ', 3E13.5 ) + 9930 FORMAT( ' WITH OUTPUT PARAMETERS:' + + / ' IERR = ', I6, + + / ' XI,ETA = ', 1P,2E13.5 ) + 9961 FORMAT( / ' ** WARNING FROM SUBROUTINE Q4_INV' + + / ' THE POINT(S) ARE OUTSIDE THE ELEMENT', + + /(' XI(',I1,') = ', 1P,2E13.5) ) + 9962 FORMAT( / ' ** WARNING FROM SUBROUTINE Q4_INV' + + / ' THE SOLUTION IS NOT UNIQUE' + + /(' XI(',I1,') = ', 1P,2E13.5) ) + 9963 FORMAT( ' CHOSE THE SOLUTION CLOSEST TO {X}' + + / ' XI,ETA = ', 1P,2E13.5, ' DELTA = ', E12.5 ) +C + END diff --git a/src/ASM/Ftn/shape2.f b/src/ASM/Ftn/shape2.f new file mode 100644 index 000000000..51907aeed --- /dev/null +++ b/src/ASM/Ftn/shape2.f @@ -0,0 +1,34 @@ + SUBROUTINE SHAPE2 (IOP,NENOD,XI,S,IPSW,IWR,IERR) + IMPLICIT NONE + INTEGER IOP,NENOD,IPSW,IWR,IERR + DOUBLE PRECISION XI(2),S(NENOD) +C + DOUBLE PRECISION UMX1,UMX2,UPX1,UPX2 +C + DOUBLE PRECISION QUART , ONE + PARAMETER ( QUART = 0.25D0, ONE = 1.0D0 ) +C + IF (IOP .EQ. 0 .AND. NENOD .EQ. 3) THEN + IERR = 0 + S(1) = XI(1) + S(2) = XI(2) + S(3) = ONE - XI(1) - XI(2) + ELSE IF (IOP .EQ. 0 .AND. NENOD .EQ. 4) THEN + IERR = 0 + UMX1 = ONE - XI(1) + UMX2 = ONE - XI(2) + UPX1 = ONE + XI(1) + UPX2 = ONE + XI(2) + S(1) = QUART * UMX1*UMX2 + S(2) = QUART * UPX1*UMX2 + S(3) = QUART * UPX1*UPX2 + S(4) = QUART * UMX1*UPX2 + ELSE + IERR = -2 + IF (IPSW .GE. 0) THEN + write(IWR,*) '*** SHAPE2:',IOP,NENOD,'not supported' + ENDIF + ENDIF +C + RETURN + END diff --git a/src/ASM/Ftn/shltri.f b/src/ASM/Ftn/shltri.f new file mode 100644 index 000000000..ceca4dc0f --- /dev/null +++ b/src/ASM/Ftn/shltri.f @@ -0,0 +1,194 @@ + SUBROUTINE SHLTRI (IOP,VX,V3,TGL,IPSW,IWR,IERR) +C +C $Id$ +C***************************** USERS' PART ***************************** +C +C A F E M ROUTINE : SHLTRI +C +C PURPOSE: +C SHLTRI : Calculate a local triad on a shell surface. +C +C METHOD: +C The local Z-axis is defined by the given vector V3. For IOP=0, the +C local X-axis is defined by the projection of the global X-axis +C onto the shell surface. If this projection is zero, the local +C Y-axis is instead defined by the projection of the global Y-axis. +C For IOP=1, the local X-axis is defined by the projection of the +C given vector VX onto the shell surface. For IOP=2, the local +C X-axis is defined by the given vector VX and the local Z-axis is +C defined by the cross product between VX and V3. +C +C INPUT ARGUMENTS: +C IOP - Option telling how to construct the triad (see above) +C VX - Vector defining the local X-axis (for IOP>0 only) +C IPSW - Print switch +C IWR - Print unit number +C +C INPUT/OUTPUT ARGUMENTS: +C V3 - Vector defining the local Z-axis +C +C OUTPUT ARGUMENTS: +C TGL - Transformation matrix +C IERR - Error flag +C +C RESTRICTIONS: +C None +C +C************************** MAINTENANCE PART *************************** +C +C INTERNAL VARIABLES: +C V1, V2 - In-plane base vectors +C +C PRINT SWITCH: +C .LT. 0 - No print +C .EQ. 0 - Print error messages only +C .GE. 6 - Print ENTER and LEAVE +C .GE. 7 - Print in addition INPUT parameters +C .GE. 8 - Print in addition OUTPUT parameters +C +C SUBROUTINES CALLED: +C DCROS3 +C DNORM3 +C +C CODING: +C 3-JUN-96, Knut Morten Okstad, NTNU +C +C REVISION: +C +C********************** END SUBROUTINE DESCRIPTION ********************* +C + IMPLICIT NONE +C +C * Global variables +C + INTEGER IOP,IPSW,IWR,IERR + DOUBLE PRECISION VX(3),V3(3),TGL(3,3) +C +C * Internal variables + INTEGER I,J +C + DOUBLE PRECISION ZERO , EPS + PARAMETER ( ZERO = 0.0D0, EPS = 1.0D-8 ) +C + DOUBLE PRECISION V1(3),V2(3) +C + DOUBLE PRECISION DNORM3 +C +C**************************** END VARIABLES **************************** +C +C Entry section +C + IF (IPSW .GE. 6) WRITE(IWR,9900) + IF (IPSW .GE. 7) THEN + WRITE(IWR,9920) IOP,V3 + IF (IOP .EQ. 1 .OR. IOP .EQ. 2) WRITE(IWR,9921) VX + ENDIF +C +C +C Program logic section +C + i = 1 + IERR = 0 + IF (IOP .NE. 2) THEN + IF (DNORM3(V3) .LE. ZERO) GOTO 7200 + ENDIF +C +C * Set up a local Cartesian basis +C (V1,V2,V3) where V3 coincide +C with the plane normal +C + IF (IOP .EQ. 1) THEN +C +C * Define V1 by the projection of VX +C onto the plane, then V2 = V3 x V1 + CALL DCROS3 (V3,VX,V2) + CALL DCROS3 (V2,V3,V1) +C + ELSE IF (IOP .EQ. 2) THEN +C + V1(1) = VX(1) + V1(2) = VX(2) + V1(3) = VX(3) +C +C * Define V3 as V1 x V3, +C then V2 = V3 x V1 + CALL DCROS3 (V1,V3,V2) + V3(1) = V2(1) + V3(2) = V2(2) + V3(3) = V2(3) + CALL DCROS3 (V3,V1,V2) + IF (DNORM3(V3) .LE. ZERO) GOTO 7200 +C + ELSE IF (ABS(V3(2)) .GT. EPS .OR. ABS(V3(3)) .GT. EPS) THEN +C +C * Define V1 by projecting the +C global X-axis onto the plane, +C then V2 = V3 x V1 +C + V1(1) = V3(2)*V3(2) + V3(3)*V3(3) + V1(2) = -V3(1)*V3(2) + V1(3) = -V3(1)*V3(3) + CALL DCROS3 (V3,V1,V2) +C + ELSE +C +C * Define V2 by projecting the +C global Y-axis onto the plane, +C then V1 = V2 x V3 + V2(1) = -V3(2)*V3(1) + V2(2) = V3(1)*V3(1) + V3(3)*V3(3) + V2(3) = -V3(2)*V3(3) + CALL DCROS3 (V2,V3,V1) +C + ENDIF +C +C * Normalize the base vectors + i = 2 + IF (DNORM3(V1) .LE. ZERO) GOTO 7200 + IF (DNORM3(V2) .LE. ZERO) GOTO 7200 +C +C * Establish the global-to-local +C transformation matrix, Tgl + DO 1000 I = 1,3 + TGL(1,I) = V1(I) + TGL(2,I) = V2(I) + TGL(3,I) = V3(I) + 1000 CONTINUE +C + GOTO 8000 +C +C Error section +C + 7200 CONTINUE + IERR = -2 + IF (IPSW .GE. 0) THEN + CALL ERRMSG ('SHLTRI',IWR,IERR) + WRITE(IWR,9920) IOP,V3 + IF (IOP .EQ. 1 .OR. IOP .EQ. 2) WRITE(IWR,9921) VX + IF (i .EQ. 2) WRITE(IWR,9922) V1,V2 + ENDIF +C +C +C Closing section +C + 8000 CONTINUE + IF (IPSW .GE. 6) WRITE(IWR,9910) + IF (IPSW .GE. 8) THEN + WRITE(IWR,9930) IERR,((TGL(i,j),j=1,3),i=1,3) + ENDIF +C + RETURN +C + 9900 FORMAT( / ' ENTERING SUBROUTINE SHLTRI' ) + 9910 FORMAT( / ' LEAVING SUBROUTINE SHLTRI' ) + 9920 FORMAT( ' WITH INPUT PARAMETERS:' + + / ' IOP = ', I6, + + / ' V3 = ', 1P,3E13.5 ) + 9921 FORMAT( ' VX = ', 1P,3E13.5 ) + 9922 FORMAT( ' V1 = ', 1P,3E13.5 + + / ' V2 = ', 3E13.5 ) + 9930 FORMAT( ' WITH OUTPUT PARAMETERS:' + + / ' IERR = ', I6, + + / ' TGL = ', 1P,3E13.5, 2(/14X,3E13.5) ) +C + END diff --git a/src/ASM/Ftn/t3_inv.f b/src/ASM/Ftn/t3_inv.f new file mode 100644 index 000000000..24dd6698c --- /dev/null +++ b/src/ASM/Ftn/t3_inv.f @@ -0,0 +1,216 @@ + SUBROUTINE T3_INV (NDIM,EPSO,EPSZ,X,TENC,XI,IPSW,IWR,IERR) +C +C $Id$ +C***************************** USERS' PART ***************************** +C +C A F E M ROUTINE : T3_INV +C +C PURPOSE: +C T3_INV : Find natural coordinates for a point within a triangle. +C +C METHOD: +C The natural (area) coordinates {XI}={A1,A2} corresponding to the +C point {X} are calculated. The element nodes are assumed to have +C the natural coordinates (1,0), (0,1) and (0,0), respectively. +C +C INPUT ARGUMENTS: +C NDIM - Number of spatial dimensions +C EPSO - Tolerence for checking if point is outside the element +C EPSZ - Tolerance for equality check with zero +C X - Cartesian coordinates of a point within the element +C TENC - Table of element nodal coordinates +C IPSW - Print switch +C IWR - Print unit number +C +C INPUT/OUTPUT ARGUMENTS: +C None +C +C OUTPUT ARGUMENTS: +C XI - Natural coordinates of the point {X} +C IERR - Error flag +C +C RESTRICTIONS: +C None +C +C************************** MAINTENANCE PART *************************** +C +C INTERNAL VARIABLES: +C AREA - Element area +C VN - Element normal vector +C V12 - Vector from node 1 to node 2 +C V13 - Vector from node 1 to node 3 +C V1X - Vector from node 1 to the point X +C +C PRINT SWITCH: +C .LT. 0 - No print +C .EQ. 0 - Print error messages only +C .GE. 7 - Print ENTER and LEAVE +C .GE. 8 - Print in addition INPUT parameters +C .GE. 9 - Print in addition OUTPUT parameters +C +C SUBROUTINES CALLED: +C DCROS3 +C DNORM3 +C +C CODING: +C 4-NOV-94, Knut Morten Okstad, NTH +C +C REVISION: +C 2-OCT-96, Knut Morten Okstad, NTNU +C The routine is included in the AFEM package. +C +C********************** END SUBROUTINE DESCRIPTION ********************* +C + IMPLICIT NONE +C +C * Global variables +C + INTEGER NDIM,IPSW,IWR,IERR + DOUBLE PRECISION EPSO,EPSZ,X(NDIM),TENC(NDIM,3),XI(2) +C +C * Internal variables +C + DOUBLE PRECISION AUX(3),A2,A3,AREA,TOL, + + V12(3),V13(3),V1X(3),VN(3) +C + DOUBLE PRECISION ZERO , ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) +C + DOUBLE PRECISION DNORM3 +C +C**************************** END VARIABLES **************************** +C +C Entry section +C + IF (IPSW .GE. 7) WRITE(IWR,9900) + IF (IPSW .GE. 8) THEN + WRITE(IWR,9920) NDIM,EPSO,EPSZ,X + CALL DPRINT (TENC,NDIM,3,'TENC ',IWR) + ENDIF +C +C +C Program logic section +C +C * Check input arguments + IERR = 0 + IF (EPSZ .LE. ZERO .OR. EPSO .LT. EPSZ) GOTO 7200 +C + IF (NDIM .EQ. 2) THEN +C +C * Calculate the element area +C + AREA = ( (TENC(1,1)-TENC(1,3))*(TENC(2,2)-TENC(2,3)) - + + (TENC(1,2)-TENC(1,3))*(TENC(2,1)-TENC(2,3)) ) + IF (AREA .LE. EPSZ) GOTO 7200 +C +C * Calculate the area coordinates +C + XI(1) = ( ( X(1) -TENC(1,3))*(TENC(2,2)-TENC(2,3)) - + + (TENC(1,2)-TENC(1,3))*( X(2) -TENC(2,3)) ) / AREA + XI(2) = ( (TENC(1,1)-TENC(1,3))*( X(2) -TENC(2,3)) - + + ( X(1) -TENC(1,3))*(TENC(2,1)-TENC(2,3)) ) / AREA +C + ELSE IF (NDIM .EQ. 3) THEN +C +C * Calculate the element normal +C vector and the element area + V12(1) = TENC(1,2) - TENC(1,1) + V12(2) = TENC(2,2) - TENC(2,1) + V12(3) = TENC(3,2) - TENC(3,1) + V13(1) = TENC(1,3) - TENC(1,1) + V13(2) = TENC(2,3) - TENC(2,1) + V13(3) = TENC(3,3) - TENC(3,1) + CALL DCROS3 (V12,V13,VN) + AREA = DNORM3(VN) + IF (AREA .LE. EPSZ) GOTO 7200 +C +C * Check the distance from point X +C to the element plane + V1X(1) = X(1) - TENC(1,1) + V1X(2) = X(2) - TENC(2,1) + V1X(3) = X(3) - TENC(3,1) + A3 = VN(1)*V1X(1) + VN(2)*V1X(2) + VN(3)*V1X(3) + IF (ABS(A3) .GT. SQRT(AREA)*EPSO) GOTO 7300 +C +C * Calculate the area coordinates + V1X(1) = V1X(1) - A3*VN(1) + V1X(2) = V1X(2) - A3*VN(2) + V1X(3) = V1X(3) - A3*VN(3) +C + CALL DCROS3 (V12,V1X,AUX) + A3 = SIGN(SQRT(AUX(1)*AUX(1)+AUX(2)*AUX(2)+AUX(3)*AUX(3)), + + AUX(1)*VN(1) +AUX(2)*VN(2) +AUX(3)*VN(3)) + CALL DCROS3 (V1X,V13,AUX) + A2 = SIGN(SQRT(AUX(1)*AUX(1)+AUX(2)*AUX(2)+AUX(3)*AUX(3)), + + AUX(1)*VN(1) +AUX(2)*VN(2) +AUX(3)*VN(3)) +C + XI(1) = (AREA-A2-A3) / AREA + XI(2) = A2 / AREA +C + ELSE + GOTO 7200 + ENDIF +C +C +C * Check that the coordinates are +C "inside" the element + TOL = ONE + EPSO + IF (XI(1) .LT. -EPSO .OR. XI(1) .GT. TOL) GOTO 6100 + IF ( XI(2) .LT. -EPSO .OR. XI(2) .GT. TOL) GOTO 6100 + IF (XI(1)+XI(2) .LT. -EPSO .OR. XI(1)+XI(2) .GT. TOL) GOTO 6100 +C + GOTO 8000 +C +C Warning section +C + 6100 CONTINUE + IERR = 1 + IF (IPSW .GE. 0) THEN + WRITE(IWR,9961) XI,X + CALL DPRINT (TENC,NDIM,3,'TENC ',IWR) + ENDIF + GOTO 8000 +C +C Error section +C + 7200 CONTINUE + IERR = -2 + GOTO 7900 + 7300 CONTINUE + IERR = -3 + IF (IPSW .GE. 0) THEN + WRITE(IWR,9973) AREA + CALL DPRINT (TENC,NDIM,3,'TENC ',IWR) + ENDIF +C + 7900 CONTINUE + IF (IPSW .GE. 0) THEN + CALL ERRMSG ('T3_INV',IWR,IERR) + ENDIF +C +C +C Closing section +C + 8000 CONTINUE + IF (IPSW .GE. 7) WRITE(IWR,9910) + IF (IPSW .GE. 9) WRITE(IWR,9930) IERR,XI +C + RETURN +C + 9900 FORMAT( / ' ENTERING SUBROUTINE T3_INV' ) + 9910 FORMAT( / ' LEAVING SUBROUTINE T3_INV' ) + 9920 FORMAT( ' WITH INPUT PARAMETERS:' + + / ' NDIM = ', I6, + + / ' EPS = ', 1P,2E13.5, + + / ' X = ', 3E13.5 ) + 9930 FORMAT( ' WITH OUTPUT PARAMETERS:' + + / ' IERR = ', I6, + + / ' XI,ETA = ', 1P,2E13.5 ) + 9961 FORMAT( / ' ** WARNING FROM SUBROUTINE T3_INV' + + / ' THE POINT IS OUTSIDE THE ELEMENT', + + / ' XI,ETA = ', 1P,2E13.5, + + / ' X = ', 3E13.5 ) + 9973 FORMAT( / ' *** THE ELEMENT IS DEGENERATED', + + / ' AREA = ', 1P,E13.5 ) +C + END diff --git a/src/ASM/Test/TestASMu2DLag.C b/src/ASM/Test/TestASMu2DLag.C index 672b2d8a8..136e67e2d 100644 --- a/src/ASM/Test/TestASMu2DLag.C +++ b/src/ASM/Test/TestASMu2DLag.C @@ -10,184 +10,164 @@ //! //============================================================================== -#include "ASMu2DLag.h" - #include "Catch2Support.h" +#include "ASMu2DLag.h" #include "MPC.h" +#include "Vec3Oper.h" +#include #include -namespace -{ -class ASMu2DLagTest : public ASMu2DLag +namespace { -public: - ASMu2DLagTest() : ASMu2DLag(2,2,'x') { ASMbase::resetNumbering(); } - - void genThreadGroups(bool separateGroup1noded = false) + std::istream& generateXMLModel(std::stringstream& str, + float Lx, float Ly, int nx, int ny, + int oneNodeElms = 0) { - this->generateThreadGroupsMultiColored(true, separateGroup1noded); - } - const ThreadGroups& getThreadGroups() const { return threadGroups; } -}; - - -void generateXMLModel(std::stringstream& str, - float Lx, float Ly, int nx, int ny, int oneNodeElms = 0) -{ - str << "\n\n"; - const double dx = Lx / nx; - const double dy = Ly / ny; - for (int j = 0; j < ny+1; ++j) - for (int i = 0; i < nx+1; ++i) - str << i*dx << " " << j*dy << " 0.0\n"; - for (int i = 0; i < oneNodeElms; ++i) - str << i+2 << " 0 0\n"; - str << "\n\n"; - for (int j = 0; j < ny; ++j) - for (int i = 0; i < nx; ++i) - str << i + j*(nx+1) << " " << i + 1 + j*(nx+1) << " " - << i + 1 + (j+1)*(nx+1) << " " << i + (j+1)*(nx+1) << '\n'; - str << "\n"; - if (oneNodeElms > 0) { - str << "\n"; + str << "\n\n"; + const double dx = Lx / nx; + const double dy = Ly / ny; + for (int j = 0; j <= ny; ++j) + for (int i = 0; i <= nx; ++i) + str << i*dx << " " << j*dy << " 0.0\n"; for (int i = 0; i < oneNodeElms; ++i) - str << i+(nx+1)*(ny+1) << '\n'; + str << i+2 << " 0 0\n"; + str << "\n\n"; + for (int j = 0; j < ny; ++j) + for (int i = 0; i < nx; ++i) + str << i + j *(nx+1) << " " << i+1 + j *(nx+1) << " " + << i + (j+1)*(nx+1) << " " << i+1 + (j+1)*(nx+1) << "\n"; str << "\n"; + if (oneNodeElms > 0) { + str << "\n"; + for (int i = 0; i < oneNodeElms; ++i) + str << (nx+1)*(ny+1) + i << "\n"; + str << "\n"; + } + str << "\n"; + return str; } - str << "\n"; -} + class ASMu2DLagTest : public ASMu2DLag + { + public: + ASMu2DLagTest(float Lx, float Ly, int nx, int ny, int oneNodeElms = 0) + : ASMu2DLag(2,2,'x') + { + ASMbase::resetNumbering(); + std::stringstream myStream; + this->read(generateXMLModel(myStream,Lx,Ly,nx,ny,oneNodeElms)); + } + + const IntMat& genThreadGroups(bool separateGroup1noded = false) + { + this->generateThreadGroupsMultiColored(false,separateGroup1noded); + return threadGroups[0]; + } + }; } TEST_CASE("TestASMu2DLag.GenerateThreadGroups3x3") { - std::stringstream str; - generateXMLModel(str, 1.0, 1.0, 3, 3); - - ASMu2DLagTest pch; - REQUIRE(pch.read(str)); + ASMu2DLagTest pch(1.0, 1.0, 3, 3); REQUIRE(pch.generateFEMTopology()); - pch.genThreadGroups(); - const ThreadGroups& groups = pch.getThreadGroups(); - REQUIRE(groups[0].size() == 4); - REQUIRE(groups[0][0].size() == 4); - REQUIRE(groups[0][1].size() == 2); - REQUIRE(groups[0][2].size() == 2); - REQUIRE(groups[0][3].size() == 1); - REQUIRE(groups[0][0] == std::vector{0, 2, 6, 8}); - REQUIRE(groups[0][1] == std::vector{1, 7}); - REQUIRE(groups[0][2] == std::vector{3, 5,}); - REQUIRE(groups[0][3] == std::vector{4}); + const IntMat& group = pch.genThreadGroups(); + + REQUIRE(group.size() == 4); + REQUIRE(group[0].size() == 4); + REQUIRE(group[1].size() == 2); + REQUIRE(group[2].size() == 2); + REQUIRE(group[3].size() == 1); + REQUIRE(group[0] == std::vector{0, 2, 6, 8}); + REQUIRE(group[1] == std::vector{1, 7}); + REQUIRE(group[2] == std::vector{3, 5,}); + REQUIRE(group[3] == std::vector{4}); } TEST_CASE("TestASMu2DLag.GenerateThreadGroups4x4") { - std::stringstream str; - generateXMLModel(str, 1.0, 1.0, 4, 4); - - ASMu2DLagTest pch; - REQUIRE(pch.read(str)); + ASMu2DLagTest pch(1.0, 1.0, 4, 4); REQUIRE(pch.generateFEMTopology()); - pch.genThreadGroups(); - const ThreadGroups& groups = pch.getThreadGroups(); - REQUIRE(groups[0].size() == 4); - REQUIRE(groups[0][0].size() == 4); - REQUIRE(groups[0][1].size() == 4); - REQUIRE(groups[0][2].size() == 4); - REQUIRE(groups[0][3].size() == 4); - REQUIRE(groups[0][0] == std::vector{0, 2, 8, 10}); - REQUIRE(groups[0][1] == std::vector{1, 3, 9, 11}); - REQUIRE(groups[0][2] == std::vector{4, 6, 12, 14}); - REQUIRE(groups[0][3] == std::vector{5, 7, 13, 15}); + const IntMat& group = pch.genThreadGroups(); + + REQUIRE(group.size() == 4); + REQUIRE(group[0].size() == 4); + REQUIRE(group[1].size() == 4); + REQUIRE(group[2].size() == 4); + REQUIRE(group[3].size() == 4); + REQUIRE(group[0] == std::vector{0, 2, 8, 10}); + REQUIRE(group[1] == std::vector{1, 3, 9, 11}); + REQUIRE(group[2] == std::vector{4, 6, 12, 14}); + REQUIRE(group[3] == std::vector{5, 7, 13, 15}); } TEST_CASE("TestASMu2DLag.GenerateThreadGroups3x3TwoPC") { - std::stringstream str; - generateXMLModel(str, 1.0, 1.0, 3, 3); - - ASMu2DLagTest pch; - REQUIRE(pch.read(str)); - pch.add2PC(1, 1, 16); - REQUIRE(pch.generateFEMTopology()); + ASMu2DLagTest pch1(1.0, 1.0, 3, 3); + pch1.add2PC(1, 1, 16); + REQUIRE(pch1.generateFEMTopology()); - str.clear(); - str.seekg(0, std::ios_base::beg); - ASMu2DLagTest pch2; - REQUIRE(pch2.read(str)); + ASMu2DLagTest pch2(1.0, 1.0, 3, 3); pch2.add2PC(16, 1, 1); REQUIRE(pch2.generateFEMTopology()); - auto checks = [](ASMu2DLagTest& p) + for (ASMu2DLagTest* pch : { &pch1, &pch2 }) { - p.genThreadGroups(); - const ThreadGroups& groups = p.getThreadGroups(); - REQUIRE(groups[0].size() == 5); - REQUIRE(groups[0][0].size() == 3); - REQUIRE(groups[0][1].size() == 2); - REQUIRE(groups[0][2].size() == 2); - REQUIRE(groups[0][3].size() == 1); - REQUIRE(groups[0][4].size() == 1); - REQUIRE(groups[0][0] == std::vector{0, 2, 6}); - REQUIRE(groups[0][1] == std::vector{1, 7}); - REQUIRE(groups[0][2] == std::vector{3, 5,}); - REQUIRE(groups[0][3] == std::vector{4}); - REQUIRE(groups[0][4] == std::vector{8}); - }; - - checks(pch); - checks(pch2); + const IntMat& group = pch->genThreadGroups(); + + REQUIRE(group.size() == 5); + REQUIRE(group[0].size() == 3); + REQUIRE(group[1].size() == 2); + REQUIRE(group[2].size() == 2); + REQUIRE(group[3].size() == 1); + REQUIRE(group[4].size() == 1); + REQUIRE(group[0] == std::vector{0, 2, 6}); + REQUIRE(group[1] == std::vector{1, 7}); + REQUIRE(group[2] == std::vector{3, 5,}); + REQUIRE(group[3] == std::vector{4}); + REQUIRE(group[4] == std::vector{8}); + } } TEST_CASE("TestASMu2DLag.GenerateThreadGroups3x3MPC") { - std::stringstream str; - generateXMLModel(str, 1.0, 1.0, 3, 3); - - ASMu2DLagTest pch; - REQUIRE(pch.read(str)); + ASMu2DLagTest pch(1.0, 1.0, 3, 3); MPC* mpc = new MPC(1, 1); mpc->addMaster(16, 1); mpc->addMaster(13, 1); pch.addMPC(mpc); REQUIRE(pch.generateFEMTopology()); - pch.genThreadGroups(); - - const ThreadGroups& groups = pch.getThreadGroups(); - REQUIRE(groups[0].size() == 4); - REQUIRE(groups[0][0].size() == 3); - REQUIRE(groups[0][1].size() == 3); - REQUIRE(groups[0][2].size() == 2); - REQUIRE(groups[0][3].size() == 1); - REQUIRE(groups[0][0] == std::vector{0, 2, 7}); - REQUIRE(groups[0][1] == std::vector{1, 6, 8}); - REQUIRE(groups[0][2] == std::vector{3, 5,}); - REQUIRE(groups[0][3] == std::vector{4}); + const IntMat& group = pch.genThreadGroups(); + + REQUIRE(group.size() == 4); + REQUIRE(group[0].size() == 3); + REQUIRE(group[1].size() == 3); + REQUIRE(group[2].size() == 2); + REQUIRE(group[3].size() == 1); + REQUIRE(group[0] == std::vector{0, 2, 7}); + REQUIRE(group[1] == std::vector{1, 6, 8}); + REQUIRE(group[2] == std::vector{3, 5,}); + REQUIRE(group[3] == std::vector{4}); } TEST_CASE("TestASMu2DLag.GenerateThreadGroups3x3OneNode") { - std::stringstream str; - generateXMLModel(str, 1.0, 1.0, 3, 3, 4); - - ASMu2DLagTest pch; - REQUIRE(pch.read(str)); + ASMu2DLagTest pch(1.0, 1.0, 3, 3, 4); REQUIRE(pch.generateFEMTopology()); - auto checks = [](ASMu2DLagTest& p, bool with1) + for (bool with1 : { false, true }) { - p.genThreadGroups(with1); - const ThreadGroups& groups = p.getThreadGroups(); + const IntMat& group = pch.genThreadGroups(with1); + const auto ref = with1 ? std::vector{ std::vector{9, 10, 11, 12}, @@ -203,32 +183,26 @@ TEST_CASE("TestASMu2DLag.GenerateThreadGroups3x3OneNode") std::vector{3, 5}, std::vector{4} }; - REQUIRE(groups[0].size() == ref.size()); + + REQUIRE(group.size() == ref.size()); for (size_t i = 0; i < ref.size(); ++i) { - REQUIRE(groups[0][i].size() == ref[i].size()); - REQUIRE(groups[0][i] == ref[i]); + REQUIRE(group[i].size() == ref[i].size()); + REQUIRE(group[i] == ref[i]); } - }; - - checks(pch, false); - checks(pch, true); + } } TEST_CASE("TestASMu2DLag.GenerateThreadGroups3x3OneNodeSPC") { - std::stringstream str; - generateXMLModel(str, 1.0, 1.0, 3, 3, 4); - - ASMu2DLagTest pch; - REQUIRE(pch.read(str)); + ASMu2DLagTest pch(1.0, 1.0, 3, 3, 4); pch.add2PC(17, 1, 1); REQUIRE(pch.generateFEMTopology()); - auto checks = [](ASMu2DLagTest& p, bool with1) + for (bool with1 : { false, true }) { - p.genThreadGroups(with1); - const ThreadGroups& groups = p.getThreadGroups(); + const IntMat& group = pch.genThreadGroups(with1); + const auto ref = with1 ? std::vector{ std::vector{10, 11, 12}, @@ -244,13 +218,82 @@ TEST_CASE("TestASMu2DLag.GenerateThreadGroups3x3OneNodeSPC") std::vector{3, 5}, std::vector{4} }; - REQUIRE(groups[0].size() == ref.size()); + + REQUIRE(group.size() == ref.size()); for (size_t i = 0; i < ref.size(); ++i) { - REQUIRE(groups[0][i].size() == ref[i].size()); - REQUIRE(groups[0][i] == ref[i]); + REQUIRE(group[i].size() == ref[i].size()); + REQUIRE(group[i] == ref[i]); } + } +} + + +TEST_CASE("TestASMu2DLag.findPoints") +{ + struct TestCase + { + const char* fname = nullptr; + int nX = 0; + int nY = 0; + size_t npt = 0; + Vec3 X; }; - checks(pch, false); - checks(pch, true); + const TestCase param = + GENERATE(TestCase{ "5x4 mesh", 5, 4, 3, Vec3(0.123,0.345,0.0) }, + TestCase{ "src/ASM/Test/refdata/Fine2Dmesh.xinp", 0, 0, 17, + Vec3(0.1,0.0,0.0) } + ); + + SECTION(param.fname) + { + std::cout <<"\nTesting "<< param.fname; + if (param.nX > 0) std::cout <<" "<< param.nX <<"x"<< param.nY; + std::cout <<" "<< param.npt << std::endl; + + ASMbase::resetNumbering(); + ASMu2DLag pch(3,3,'x'); + if (param.nX > 0 && param.nY > 0) + { + // Generate a simple mesh on the bi-unit square + std::stringstream str; + REQUIRE(pch.read(generateXMLModel(str,1.0,1.0,param.nX,param.nY))); + } + else + { + // Read mesh from file + std::ifstream fs(param.fname); + REQUIRE(pch.read(fs)); + } + REQUIRE(pch.generateFEMTopology()); + + // Pick some random elements and search for their center points + const size_t nel = pch.getNoElms(); + const size_t npt = param.npt; + std::vector elements; elements.reserve(npt); + std::vector points; points.reserve(npt); + for (size_t iel = 1+nel/(2*npt); iel <= nel; iel += nel/npt) + { + elements.push_back(iel); + points.emplace_back(pch.getElementCenter(iel)); + std::cout <<"Element "<< iel <<": "<< points.back() << std::endl; + } + points.emplace_back(param.X); + + std::vector elms; + REQUIRE(pch.findPoints(points,elms)); + + size_t ipt = 0; + std::cout <<"Search results:"; + for (const ASMbase::PointParams& elm : elms) + { + std::cout <<"\nPoint "<< ++ipt <<": iel = "<< elm.iel + <<" xi,eta = "<< elm.u[0] <<","<< elm.u[1] + <<" dist = "<< elm.dist; + if (ipt-1 < elements.size()) + REQUIRE(elements[ipt-1] == elm.iel); + REQUIRE(elm.dist <= 1.0e-12); + } + std::cout << std::endl; + } } diff --git a/src/ASM/Test/refdata/Fine2Dmesh.xinp b/src/ASM/Test/refdata/Fine2Dmesh.xinp new file mode 100644 index 000000000..8f5c832b2 --- /dev/null +++ b/src/ASM/Test/refdata/Fine2Dmesh.xinp @@ -0,0 +1,1861 @@ + + + 0.100488 0 -0.00270879 + 0.100458 0 -0.00246075 + 0.100423 0 -0.00226954 + 0.100375 0 -0.00207408 + 0.100287 0 -0.00179642 + 0.100175 0 -0.00152714 + 0.100041 0 -0.00126816 + 0.0999368 0 -0.00109605 + 0.0998266 0 -0.000936094 + 0.0996721 0 -0.00073968 + 0.0994749 0 -0.000525124 + 0.0992603 0 -0.000327851 + 0.0990639 0 -0.000173401 + 0.0989039 0 -6.31507e-05 + 0.0987318 0 4.12858e-05 + 0.0984729 0 0.000174987 + 0.0982036 0 0.000286524 + 0.0979259 0 0.000375111 + 0.0977304 0 0.000422931 + 0.0975393 0 0.000458118 + 0.0972912 0 0.000487782 + 0.100987 0 -0.00267906 + 0.100954 0 -0.00239547 + 0.100911 0 -0.00216132 + 0.100852 0 -0.00192463 + 0.100754 0 -0.00161889 + 0.100631 0 -0.00132213 + 0.100484 0 -0.00103632 + 0.100359 0 -0.000827405 + 0.100223 0 -0.00063166 + 0.100046 0 -0.000407722 + 0.0998284 0 -0.000171571 + 0.0995923 0 4.61626e-05 + 0.0993683 0 0.00022334 + 0.0991726 0 0.000358553 + 0.0989637 0 0.000484495 + 0.0986779 0 0.000631049 + 0.0983811 0 0.00075397 + 0.0980754 0 0.000852446 + 0.0978386 0 0.000911088 + 0.0976046 0 0.000953922 + 0.0973209 0 0.000986979 + 0.101486 0 -0.00264932 + 0.10145 0 -0.00233019 + 0.101399 0 -0.00205311 + 0.10133 0 -0.00177517 + 0.101221 0 -0.00144137 + 0.101087 0 -0.00111712 + 0.100928 0 -0.000804473 + 0.10078 0 -0.000558757 + 0.10062 0 -0.000327226 + 0.10042 0 -7.57639e-05 + 0.100182 0 0.000181982 + 0.0999242 0 0.000420176 + 0.0996728 0 0.000620082 + 0.0994412 0 0.000780257 + 0.0991955 0 0.000927705 + 0.0988829 0 0.00108711 + 0.0985586 0 0.00122142 + 0.0982248 0 0.00132978 + 0.0979468 0 0.00139924 + 0.0976698 0 0.00144973 + 0.0973507 0 0.00148618 + 0.101985 0 -0.00261958 + 0.101946 0 -0.00226491 + 0.101887 0 -0.0019449 + 0.101807 0 -0.00162571 + 0.101689 0 -0.00126384 + 0.101543 0 -0.000912119 + 0.101371 0 -0.00057263 + 0.101202 0 -0.000290108 + 0.101017 0 -2.27914e-05 + 0.100794 0 0.000256194 + 0.100536 0 0.000535535 + 0.100256 0 0.00079419 + 0.0999772 0 0.00101682 + 0.0997099 0 0.00120196 + 0.0994274 0 0.00137091 + 0.0990879 0 0.00154317 + 0.0987362 0 0.00168886 + 0.0983743 0 0.00180711 + 0.0980551 0 0.0018874 + 0.0977351 0 0.00194553 + 0.0973804 0 0.00198537 + 0.102485 0 -0.00258984 + 0.102441 0 -0.00219963 + 0.102376 0 -0.00183668 + 0.102284 0 -0.00147625 + 0.102156 0 -0.00108632 + 0.101999 0 -0.000707113 + 0.101814 0 -0.000340786 + 0.101624 0 -2.14593e-05 + 0.101414 0 0.000281643 + 0.101168 0 0.000588152 + 0.100889 0 0.000889088 + 0.100588 0 0.0011682 + 0.100282 0 0.00141356 + 0.0999785 0 0.00162366 + 0.0996592 0 0.00181412 + 0.0992929 0 0.00199924 + 0.0989137 0 0.00215631 + 0.0985238 0 0.00228445 + 0.0981633 0 0.00237556 + 0.0978004 0 0.00244134 + 0.0974102 0 0.00248457 + 0.102984 0 -0.0025601 + 0.102937 0 -0.00213435 + 0.102864 0 -0.00172847 + 0.102762 0 -0.00132679 + 0.102624 0 -0.000908797 + 0.102455 0 -0.000502107 + 0.102257 0 -0.000108943 + 0.102045 0 0.000247189 + 0.10181 0 0.000586077 + 0.101542 0 0.00092011 + 0.101243 0 0.00124264 + 0.10092 0 0.00154222 + 0.100586 0 0.00181031 + 0.100247 0 0.00204537 + 0.0998911 0 0.00225733 + 0.0994979 0 0.0024553 + 0.0990912 0 0.00262376 + 0.0986732 0 0.00276178 + 0.0982715 0 0.00286371 + 0.0978657 0 0.00293714 + 0.0974399 0 0.00298377 + 0.103483 0 -0.00253036 + 0.103433 0 -0.00206908 + 0.103352 0 -0.00162026 + 0.103239 0 -0.00117733 + 0.103091 0 -0.000731274 + 0.102911 0 -0.000297101 + 0.102701 0 0.000122901 + 0.102467 0 0.000515838 + 0.102207 0 0.000890511 + 0.101916 0 0.00125207 + 0.101596 0 0.00159619 + 0.101252 0 0.00191623 + 0.100891 0 0.00220705 + 0.100516 0 0.00246707 + 0.100123 0 0.00270054 + 0.0997029 0 0.00291136 + 0.0992687 0 0.0030912 + 0.0988227 0 0.00323912 + 0.0983797 0 0.00335187 + 0.0979309 0 0.00343295 + 0.0974696 0 0.00348297 + 0.097 0 0.0035 + 0.097 0 0.003 + 0.097 0 0.0025 + 0.097 0 0.002 + 0.097 0 0.0015 + 0.097 0 0.001 + 0.097 0 0.0005 + 0.0972615 0 -1.14159e-05 + 0.097474 0 -3.76868e-05 + 0.0976222 0 -6.52253e-05 + 0.0977765 0 -0.000102223 + 0.0980261 0 -0.000180923 + 0.0982679 0 -0.000281076 + 0.0985 0 -0.000401924 + 0.0986352 0 -0.000484854 + 0.0987595 0 -0.000570142 + 0.0989284 0 -0.000701865 + 0.0991213 0 -0.000878677 + 0.0992981 0 -0.00107164 + 0.0994299 0 -0.00124053 + 0.0995151 0 -0.0013647 + 0.0995981 0 -0.0015 + 0.0997189 0 -0.00173214 + 0.0998191 0 -0.00197394 + 0.0998978 0 -0.00222354 + 0.0999348 0 -0.00237775 + 0.0999623 0 -0.00252603 + 0.0999886 0 -0.00273853 + 0.103982 0 -0.00250063 + 0.103929 0 -0.0020038 + 0.10384 0 -0.00151204 + 0.103716 0 -0.00102787 + 0.103559 0 -0.00055375 + 0.103367 0 -9.20946e-05 + 0.103144 0 0.000354744 + 0.102889 0 0.000784487 + 0.102604 0 0.00119495 + 0.10229 0 0.00158403 + 0.10195 0 0.00194975 + 0.101584 0 0.00229025 + 0.101195 0 0.00260379 + 0.100784 0 0.00288877 + 0.100355 0 0.00314375 + 0.0999079 0 0.00336742 + 0.0994462 0 0.00355865 + 0.0989721 0 0.00371645 + 0.098488 0 0.00384003 + 0.0979962 0 0.00392875 + 0.0974994 0 0.00398216 + 0.1035 0 -0.003 + 0.103 0 -0.003 + 0.1025 0 -0.003 + 0.102 0 -0.003 + 0.1015 0 -0.003 + 0.101 0 -0.003 + 0.1005 0 -0.003 + 0.097 0 0 + 0.097 0 0.004 + 0.1 0 -0.003 + 0.104 0 -0.003 + 0.12665 0 0.0289971 + 0.125661 0 0.0289445 + 0.124446 0 0.0288412 + 0.123023 0 0.0287277 + 0.12148 0 0.0285904 + 0.119865 0 0.0283934 + 0.11806 0 0.0283224 + 0.116262 0 0.0282936 + 0.114485 0 0.0282789 + 0.112723 0 0.0282689 + 0.110969 0 0.0282602 + 0.10922 0 0.0282521 + 0.107472 0 0.0282442 + 0.105726 0 0.0282367 + 0.10398 0 0.0282296 + 0.102235 0 0.0282231 + 0.10049 0 0.0282171 + 0.0987449 0 0.0282118 + 0.0987405 0 0.026423 + 0.0987342 0 0.0246333 + 0.0987252 0 0.0228425 + 0.0987135 0 0.0210504 + 0.0986995 0 0.0192565 + 0.0986837 0 0.0174593 + 0.0986664 0 0.0156552 + 0.0986454 0 0.0138387 + 0.0986115 0 0.0120046 + 0.09853 0 0.0101753 + 0.098273 0 0.00856521 + 0.0980445 0 0.00722559 + 0.0978595 0 0.0061374 + 0.0977098 0 0.00527673 + 0.0975947 0 0.00457822 + 0.0981704 0 0.00452173 + 0.0987063 0 0.00439854 + 0.0991867 0 0.00425457 + 0.0996818 0 0.00406872 + 0.100156 0 0.00385058 + 0.100595 0 0.003631 + 0.101063 0 0.00337721 + 0.101526 0 0.00307077 + 0.101952 0 0.00272058 + 0.102352 0 0.00235359 + 0.102687 0 0.00194511 + 0.102997 0 0.00155958 + 0.103329 0 0.00114126 + 0.10363 0 0.000666665 + 0.103862 0 0.000173918 + 0.104053 0 -0.000317079 + 0.104215 0 -0.000798186 + 0.104378 0 -0.00128132 + 0.104517 0 -0.00181102 + 0.104579 0 -0.0023927 + 0.105284 0 -0.00225575 + 0.106168 0 -0.0020612 + 0.107303 0 -0.00180852 + 0.108726 0 -0.00169339 + 0.110311 0 -0.00155966 + 0.11205 0 -0.00146469 + 0.113841 0 -0.00140362 + 0.115644 0 -0.00136419 + 0.117449 0 -0.00133748 + 0.11925 0 -0.0013181 + 0.121049 0 -0.00130295 + 0.122844 0 -0.00129034 + 0.124636 0 -0.00127939 + 0.126426 0 -0.00126962 + 0.128213 0 -0.00126058 + 0.12822 0 0.00047981 + 0.128226 0 0.00222217 + 0.128231 0 0.00396681 + 0.128236 0 0.00571375 + 0.128241 0 0.00746276 + 0.128246 0 0.0092137 + 0.128252 0 0.010967 + 0.12826 0 0.0127245 + 0.12827 0 0.0144899 + 0.128286 0 0.01627 + 0.128319 0 0.0180724 + 0.128397 0 0.0198807 + 0.128602 0 0.0214944 + 0.128746 0 0.023027 + 0.128852 0 0.0244135 + 0.128902 0 0.025545 + 0.12891 0 0.0265387 + 0.128815 0 0.0274124 + 0.128556 0 0.0281426 + 0.128133 0 0.0286724 + 0.127481 0 0.0289285 + 0.126466 0 0.0281245 + 0.12554 0 0.0279441 + 0.124385 0 0.027695 + 0.123013 0 0.0274789 + 0.121569 0 0.0272373 + 0.120205 0 0.026656 + 0.118144 0 0.0265956 + 0.116277 0 0.0265654 + 0.114476 0 0.0265459 + 0.112705 0 0.0265298 + 0.110947 0 0.0265145 + 0.109197 0 0.0264994 + 0.10745 0 0.0264846 + 0.105706 0 0.0264703 + 0.103963 0 0.0264568 + 0.102222 0 0.0264444 + 0.100481 0 0.0264331 + 0.100469 0 0.0246475 + 0.100452 0 0.0228596 + 0.10043 0 0.0210693 + 0.100405 0 0.0192759 + 0.100377 0 0.0174767 + 0.100348 0 0.0156658 + 0.100318 0 0.0138297 + 0.100285 0 0.0119321 + 0.100249 0 0.00985007 + 0.0995275 0 0.00824778 + 0.0990525 0 0.00698662 + 0.0986878 0 0.00595154 + 0.0983874 0 0.00518214 + 0.0989839 0 0.00497915 + 0.0993911 0 0.00483662 + 0.0999436 0 0.00461402 + 0.100445 0 0.0043389 + 0.100815 0 0.00415674 + 0.101353 0 0.00392321 + 0.101903 0 0.00358705 + 0.102354 0 0.00317382 + 0.102826 0 0.00280539 + 0.103114 0 0.00228998 + 0.103373 0 0.00196416 + 0.103806 0 0.00156255 + 0.104187 0 0.00100622 + 0.1044 0 0.000448633 + 0.10458 0 -7.92453e-05 + 0.104715 0 -0.000555746 + 0.104944 0 -0.000994438 + 0.105183 0 -0.00155905 + 0.10599 0 -0.00116758 + 0.10711 0 -0.000464572 + 0.108683 0 -0.000390918 + 0.110195 0 -6.73027e-05 + 0.111985 0 0.000115755 + 0.113812 0 0.000223629 + 0.115636 0 0.000292083 + 0.117452 0 0.000339106 + 0.119261 0 0.000374 + 0.121063 0 0.000401768 + 0.122859 0 0.000425108 + 0.12465 0 0.000445418 + 0.126436 0 0.000463366 + 0.126447 0 0.00219949 + 0.126457 0 0.00393984 + 0.126467 0 0.00568469 + 0.126476 0 0.00743383 + 0.126487 0 0.00918698 + 0.126499 0 0.0109447 + 0.126512 0 0.0127094 + 0.126529 0 0.0144874 + 0.126552 0 0.0162936 + 0.126591 0 0.0181674 + 0.12667 0 0.0202362 + 0.127266 0 0.0216007 + 0.127533 0 0.0230384 + 0.127761 0 0.0243603 + 0.127849 0 0.0252196 + 0.127923 0 0.0262034 + 0.127939 0 0.0270476 + 0.12785 0 0.0277287 + 0.127732 0 0.0281975 + 0.127213 0 0.0281954 + 0.123089 0 0.0262564 + 0.124556 0 0.0265196 + 0.127174 0 0.0275357 + 0.102204 0 0.0246633 + 0.103941 0 0.0246809 + 0.105681 0 0.0247 + 0.107423 0 0.0247201 + 0.109168 0 0.0247408 + 0.110917 0 0.0247613 + 0.112674 0 0.024781 + 0.114444 0 0.0247995 + 0.116233 0 0.0248165 + 0.118048 0 0.02483 + 0.119872 0 0.0248264 + 0.121525 0 0.0247355 + 0.123211 0 0.0248122 + 0.125098 0 0.025117 + 0.102181 0 0.022879 + 0.103914 0 0.0229008 + 0.10565 0 0.0229245 + 0.10739 0 0.0229495 + 0.109133 0 0.0229749 + 0.110882 0 0.0229998 + 0.112637 0 0.0230232 + 0.114401 0 0.0230443 + 0.116175 0 0.0230624 + 0.117956 0 0.0230766 + 0.119726 0 0.0230846 + 0.121441 0 0.0230872 + 0.123143 0 0.0231488 + 0.124838 0 0.0232475 + 0.102152 0 0.0210915 + 0.10388 0 0.0211168 + 0.105613 0 0.0211445 + 0.107351 0 0.0211736 + 0.109094 0 0.0212029 + 0.110842 0 0.0212316 + 0.112596 0 0.0212586 + 0.114356 0 0.0212833 + 0.116122 0 0.0213056 + 0.117888 0 0.021327 + 0.119645 0 0.0213511 + 0.121378 0 0.0213876 + 0.123089 0 0.0214589 + 0.124756 0 0.0215609 + 0.102119 0 0.0193003 + 0.103841 0 0.019329 + 0.105571 0 0.0193606 + 0.107308 0 0.0193935 + 0.10905 0 0.0194266 + 0.110799 0 0.0194588 + 0.112552 0 0.0194894 + 0.114311 0 0.0195182 + 0.116074 0 0.0195457 + 0.117837 0 0.0195745 + 0.119595 0 0.0196099 + 0.121342 0 0.0196621 + 0.123083 0 0.0197495 + 0.124833 0 0.0199048 + 0.102082 0 0.0175034 + 0.103799 0 0.0175367 + 0.105525 0 0.0175731 + 0.107261 0 0.0176104 + 0.109003 0 0.0176472 + 0.110752 0 0.0176829 + 0.112507 0 0.0177171 + 0.114267 0 0.01775 + 0.11603 0 0.0177824 + 0.117794 0 0.0178167 + 0.119556 0 0.0178571 + 0.121313 0 0.0179098 + 0.123066 0 0.0179822 + 0.124823 0 0.0180772 + 0.102045 0 0.0156972 + 0.103755 0 0.0157387 + 0.105478 0 0.0157823 + 0.107211 0 0.015825 + 0.108953 0 0.0158659 + 0.110702 0 0.0159051 + 0.112459 0 0.0159429 + 0.114221 0 0.0159796 + 0.115987 0 0.0160162 + 0.117754 0 0.0160542 + 0.119521 0 0.0160961 + 0.121284 0 0.0161445 + 0.123043 0 0.0162001 + 0.124799 0 0.0162578 + 0.10201 0 0.0138763 + 0.103713 0 0.0139345 + 0.105429 0 0.0139892 + 0.107158 0 0.0140387 + 0.108899 0 0.014084 + 0.110649 0 0.0141265 + 0.112408 0 0.0141675 + 0.114173 0 0.0142077 + 0.115943 0 0.0142479 + 0.117715 0 0.0142888 + 0.119486 0 0.0143314 + 0.121255 0 0.0143759 + 0.123018 0 0.0144212 + 0.124776 0 0.0144624 + 0.101982 0 0.012038 + 0.103675 0 0.0121291 + 0.10538 0 0.0121979 + 0.107103 0 0.0122539 + 0.108841 0 0.0123028 + 0.110591 0 0.012348 + 0.112353 0 0.0123917 + 0.114122 0 0.0124351 + 0.115897 0 0.0124786 + 0.117674 0 0.0125221 + 0.119451 0 0.0125654 + 0.121225 0 0.0126081 + 0.122994 0 0.0126486 + 0.124756 0 0.0126842 + 0.101979 0 0.0102086 + 0.10364 0 0.0103413 + 0.105327 0 0.0104162 + 0.107042 0 0.0104743 + 0.108776 0 0.0105248 + 0.110528 0 0.0105711 + 0.112293 0 0.0106168 + 0.114068 0 0.010663 + 0.115849 0 0.0107096 + 0.117633 0 0.0107556 + 0.119417 0 0.0108002 + 0.121197 0 0.0108426 + 0.122971 0 0.010882 + 0.124737 0 0.0109167 + 0.102064 0 0.00860037 + 0.103598 0 0.00860859 + 0.105262 0 0.00864992 + 0.106971 0 0.00870297 + 0.108703 0 0.0087527 + 0.110455 0 0.00879785 + 0.112226 0 0.00884427 + 0.114009 0 0.00889317 + 0.1158 0 0.00894289 + 0.117592 0 0.00899148 + 0.119384 0 0.0090377 + 0.12117 0 0.00908094 + 0.122949 0 0.00912069 + 0.124721 0 0.00915629 + 0.101756 0 0.0069012 + 0.103443 0 0.00684337 + 0.105163 0 0.00687461 + 0.106886 0 0.00693964 + 0.108616 0 0.00699155 + 0.110369 0 0.00703038 + 0.11215 0 0.00707616 + 0.113947 0 0.0071283 + 0.11575 0 0.00718182 + 0.117553 0 0.00723325 + 0.119352 0 0.00728114 + 0.121145 0 0.0073252 + 0.122929 0 0.00736544 + 0.124705 0 0.00740179 + 0.101965 0 0.00546758 + 0.103269 0 0.0049944 + 0.105072 0 0.00507335 + 0.106804 0 0.00519968 + 0.108513 0 0.00525533 + 0.11026 0 0.00526811 + 0.112064 0 0.00531278 + 0.113883 0 0.00537142 + 0.115702 0 0.0054306 + 0.117516 0 0.00548516 + 0.119322 0 0.00553414 + 0.12112 0 0.00557806 + 0.12291 0 0.00561762 + 0.12469 0 0.00565332 + 0.103493 0 0.00341329 + 0.105061 0 0.00323793 + 0.106754 0 0.00354626 + 0.108377 0 0.00357906 + 0.110101 0 0.00349043 + 0.111968 0 0.00354996 + 0.113825 0 0.00362628 + 0.115662 0 0.00369438 + 0.117485 0 0.00375168 + 0.119296 0 0.00380007 + 0.121098 0 0.00384182 + 0.122891 0 0.00387858 + 0.124676 0 0.00391144 + 0.105658 0 0.00203031 + 0.106778 0 0.00218835 + 0.108146 0 0.00204454 + 0.109806 0 0.0015831 + 0.111887 0 0.00178606 + 0.11379 0 0.00190254 + 0.115638 0 0.00198019 + 0.117463 0 0.0020369 + 0.119275 0 0.00208118 + 0.121078 0 0.00211761 + 0.122873 0 0.00214883 + 0.124661 0 0.00217629 + 0.101643 0 0.00457787 + 0.100793 0 0.00600991 + 0.100265 0 0.00521521 + 0.101209 0 0.00540471 + 0.100951 0 0.00867966 + 0.09996 0 0.00652133 + 0.100549 0 0.00758749 + 0.102744 0 0.00359039 + 0.10434 0 0.0021509 + 0.107629 0 0.000848334 + 0.105156 0 0.000123113 + 0.10595 0 0.000225313 + 0.106559 0 0.00115345 + 0.126361 0 0.0231743 + 0.10239 0 0.00418743 + 0.105882 0 0.00136163 + 0.105124 0 -0.000336288 + 0.105502 0 -0.000568249 + 0.126466 0 0.0273734 + 0.12566 0 0.0270188 + 0.108706 0 0.000680171 + 0.126995 0 0.0260161 + 0.127169 0 0.0268546 + 0.126081 0 0.0262427 + 0.126572 0 0.0268236 + 0.121597 0 0.0262213 + 0.100875 0 0.00473868 + 0.0994576 0 0.00550077 + 0.10358 0 0.00245458 + 0.10543 0 0.000910276 + 0.104975 0 0.000708883 + 0.104918 0 0.00136132 + 0.126263 0 0.021643 + 0.126813 0 0.0247775 + 0.097 0 0.0282067 + 0.097 0 0.0264134 + 0.097 0 0.0246201 + 0.097 0 0.0228269 + 0.097 0 0.0210336 + 0.097 0 0.0192403 + 0.097 0 0.017447 + 0.097 0 0.0156537 + 0.097 0 0.0138604 + 0.097 0 0.0120672 + 0.097 0 0.0102772 + 0.097 0 0.00860989 + 0.097 0 0.00721229 + 0.097 0 0.00609581 + 0.097 0 0.00520954 + 0.097 0 0.0045325 + 0.126014 0 0.03 + 0.124727 0 0.03 + 0.123168 0 0.03 + 0.121476 0 0.03 + 0.119728 0 0.03 + 0.117979 0 0.03 + 0.116231 0 0.03 + 0.114483 0 0.03 + 0.112735 0 0.03 + 0.110986 0 0.03 + 0.109238 0 0.03 + 0.10749 0 0.03 + 0.105741 0 0.03 + 0.103993 0 0.03 + 0.102245 0 0.03 + 0.100497 0 0.03 + 0.0987483 0 0.03 + 0.128207 0 -0.003 + 0.126413 0 -0.003 + 0.12462 0 -0.003 + 0.122827 0 -0.003 + 0.121034 0 -0.003 + 0.11924 0 -0.003 + 0.117447 0 -0.003 + 0.115654 0 -0.003 + 0.11386 0 -0.003 + 0.112067 0 -0.003 + 0.110277 0 -0.003 + 0.10861 0 -0.003 + 0.107212 0 -0.003 + 0.106096 0 -0.003 + 0.10521 0 -0.003 + 0.104533 0 -0.003 + 0.129853 0 0.027927 + 0.129427 0 0.0287634 + 0.128763 0 0.029427 + 0.127927 0 0.0298532 + 0.13 0 -0.00125171 + 0.13 0 0.000496575 + 0.13 0 0.00224486 + 0.13 0 0.00399315 + 0.13 0 0.00574144 + 0.13 0 0.00748972 + 0.13 0 0.00923801 + 0.13 0 0.0109863 + 0.13 0 0.0127346 + 0.13 0 0.0144829 + 0.13 0 0.0162312 + 0.13 0 0.0179794 + 0.13 0 0.0197277 + 0.13 0 0.021476 + 0.13 0 0.0231677 + 0.13 0 0.0247271 + 0.13 0 0.0260137 + 0.097 0 0.03 + 0.127 0 0.03 + 0.13 0 -0.003 + 0.13 0 0.027 + 0.0884378 0 0.0281585 + 0.0901539 0 0.0281953 + 0.0918678 0 0.0282051 + 0.09358 0 0.0282077 + 0.0952908 0 0.0282077 + 0.0884483 0 0.0264379 + 0.0901657 0 0.0264177 + 0.0918785 0 0.0264174 + 0.0935878 0 0.0264178 + 0.0952946 0 0.0264163 + 0.0867357 0 0.0245771 + 0.0884638 0 0.024603 + 0.0901823 0 0.0246203 + 0.0918929 0 0.0246288 + 0.093598 0 0.02463 + 0.0952998 0 0.0246263 + 0.0884889 0 0.0227769 + 0.090207 0 0.0228318 + 0.0919128 0 0.0228475 + 0.0936115 0 0.0228471 + 0.0953065 0 0.0228388 + 0.088532 0 0.021085 + 0.0902441 0 0.0210823 + 0.0919399 0 0.0210821 + 0.0936287 0 0.0210722 + 0.0953148 0 0.0210547 + 0.0868662 0 0.0192233 + 0.0886073 0 0.0193021 + 0.0902975 0 0.0193303 + 0.091974 0 0.0193264 + 0.0936487 0 0.0193049 + 0.0953241 0 0.0192744 + 0.0870697 0 0.0173217 + 0.0887337 0 0.0175696 + 0.0903646 0 0.0176103 + 0.09201 0 0.0175883 + 0.0936681 0 0.0175464 + 0.0953328 0 0.0174977 + 0.0876791 0 0.0161301 + 0.0888933 0 0.0160442 + 0.0904174 0 0.0159531 + 0.0920335 0 0.0158701 + 0.0936809 0 0.0157944 + 0.095339 0 0.0157232 + 0.0870745 0 0.0150243 + 0.0887432 0 0.0145238 + 0.0903782 0 0.0142874 + 0.0920259 0 0.0141444 + 0.0936831 0 0.0140381 + 0.0953424 0 0.0139468 + 0.086876 0 0.0130795 + 0.0886271 0 0.012739 + 0.0903264 0 0.012528 + 0.092009 0 0.0123817 + 0.0936832 0 0.0122666 + 0.0953472 0 0.012166 + 0.0868032 0 0.0109182 + 0.0885635 0 0.0108242 + 0.0902921 0 0.0107037 + 0.0920009 0 0.0105875 + 0.0936922 0 0.0104828 + 0.0953618 0 0.0103857 + 0.0885347 0 0.00893657 + 0.0902799 0 0.00887591 + 0.0920082 0 0.00878814 + 0.0937173 0 0.00870094 + 0.0954037 0 0.0086223 + 0.0885292 0 0.0071904 + 0.090286 0 0.0070862 + 0.0920265 0 0.00700692 + 0.0937526 0 0.00692854 + 0.0955286 0 0.0068007 + 0.0867698 0 0.00540108 + 0.0885411 0 0.00534604 + 0.0903062 0 0.00529404 + 0.0920444 0 0.00525206 + 0.0937195 0 0.0052281 + 0.0952284 0 0.00522512 + 0.0885675 0 0.00351718 + 0.090348 0 0.00351661 + 0.0921082 0 0.00350542 + 0.0938315 0 0.00352788 + 0.0954514 0 0.00369025 + 0.0886041 0 0.00181618 + 0.0904041 0 0.00176717 + 0.0921952 0 0.00174188 + 0.0960787 0 0.00147056 + 0.0965341 0 0.00184952 + 0.0963531 0 0.00243252 + 0.0961985 0 0.00519451 + 0.0963132 0 0.00432189 + 0.0962424 0 0.00603034 + 0.0964519 0 0.000612568 + 0.095818 0 0.000807326 + 0.0965991 0 0.00387286 + 0.0965576 0 0.00114515 + 0.0966972 0 0.00149822 + 0.0867196 0 0.0280008 + 0.0867259 0 0.0265721 + 0.0867531 0 0.0225877 + 0.0867878 0 0.0211785 + 0.0867758 0 0.00886077 + 0.0867687 0 0.00740852 + 0.0867846 0 0.00340786 + 0.0868026 0 0.00198783 + 0.0940323 0 0.00169997 + 0.0953878 0 0.0020879 + 0.0950197 0 0.0011488 + 0.0964807 0 0.00329913 + 0.103638 0 -0.0046441 + 0.105189 0 -0.00481537 + 0.106779 0 -0.00449506 + 0.108593 0 -0.00461211 + 0.11034 0 -0.0046468 + 0.112103 0 -0.00465749 + 0.113872 0 -0.00466331 + 0.115648 0 -0.00467149 + 0.117432 0 -0.00468258 + 0.119225 0 -0.00469506 + 0.121024 0 -0.00470759 + 0.122825 0 -0.00471919 + 0.124627 0 -0.0047284 + 0.126426 0 -0.00473243 + 0.128215 0 -0.00472564 + 0.103509 0 -0.00622966 + 0.105205 0 -0.00633575 + 0.106897 0 -0.00628743 + 0.108644 0 -0.0063066 + 0.110387 0 -0.00631763 + 0.112132 0 -0.00631988 + 0.113877 0 -0.00632435 + 0.115633 0 -0.00634014 + 0.117409 0 -0.00636385 + 0.119204 0 -0.00639016 + 0.121011 0 -0.00641622 + 0.122823 0 -0.00644088 + 0.124637 0 -0.00646212 + 0.126439 0 -0.00646766 + 0.128223 0 -0.00645166 + 0.101737 0 -0.00783441 + 0.103493 0 -0.00794336 + 0.105231 0 -0.00800666 + 0.106963 0 -0.00801005 + 0.108699 0 -0.00800812 + 0.110434 0 -0.0079968 + 0.11216 0 -0.00797995 + 0.113872 0 -0.00797409 + 0.115598 0 -0.00800096 + 0.117367 0 -0.00804261 + 0.119172 0 -0.00808571 + 0.120993 0 -0.0081267 + 0.122818 0 -0.00816657 + 0.124644 0 -0.0081964 + 0.126447 0 -0.00819828 + 0.128224 0 -0.00817363 + 0.10176 0 -0.00962154 + 0.103504 0 -0.00969011 + 0.105269 0 -0.00972834 + 0.107028 0 -0.00973266 + 0.108756 0 -0.00971642 + 0.110493 0 -0.00968036 + 0.112202 0 -0.00962863 + 0.113854 0 -0.00959102 + 0.115521 0 -0.00964705 + 0.117288 0 -0.00972021 + 0.119125 0 -0.00978428 + 0.120974 0 -0.00983971 + 0.1228 0 -0.00988831 + 0.124636 0 -0.00991727 + 0.126447 0 -0.00991438 + 0.128213 0 -0.00988758 + 0.101811 0 -0.0114113 + 0.103507 0 -0.0114518 + 0.105321 0 -0.0114725 + 0.10713 0 -0.0114691 + 0.108808 0 -0.0114411 + 0.110579 0 -0.0113781 + 0.112302 0 -0.0112626 + 0.11382 0 -0.0111141 + 0.115346 0 -0.0112762 + 0.11714 0 -0.0114078 + 0.119066 0 -0.0114933 + 0.120972 0 -0.0115532 + 0.122742 0 -0.0115961 + 0.124611 0 -0.0116185 + 0.126457 0 -0.0116145 + 0.12817 0 -0.0115938 + 0.105384 0 -0.0132316 + 0.110716 0 -0.0131268 + 0.112608 0 -0.012929 + 0.113777 0 -0.0123264 + 0.114905 0 -0.0129363 + 0.116859 0 -0.0131428 + 0.119028 0 -0.0132271 + 0.124577 0 -0.0133094 + 0.103224 0 -0.00397825 + 0.105138 0 -0.00378865 + 0.106004 0 -0.00376124 + 0.101528 0 -0.00370832 + 0.102597 0 -0.00382751 + 0.101039 0 -0.00341169 + 0.10082 0 -0.0041226 + 0.103439 0 -0.00346991 + 0.10059 0 -0.0035272 + 0.101984 0 -0.0132038 + 0.103401 0 -0.0132218 + 0.107366 0 -0.0132251 + 0.108768 0 -0.0131983 + 0.121096 0 -0.0132702 + 0.122558 0 -0.0132939 + 0.12658 0 -0.0133066 + 0.128006 0 -0.0132964 + 0.104149 0 -0.00358053 + 0.102156 0 -0.00338396 + 0.102884 0 -0.00332436 + 0.101162 0 -0.00497468 + 0.102126 0 -0.00462898 + 0.101707 0 -0.00599016 + 0.085 0 0.0272727 + 0.085 0 0.0245455 + 0.085 0 0.0218182 + 0.085 0 0.0190909 + 0.085 0 0.0163636 + 0.085 0 0.0136364 + 0.085 0 0.0109091 + 0.085 0 0.00818182 + 0.085 0 0.00545455 + 0.085 0 0.00272727 + 0.0868219 0 0 + 0.0886438 0 0 + 0.0904656 0 0 + 0.0922304 0 0 + 0.0937169 0 0 + 0.0948407 0 0 + 0.0957179 0 0 + 0.0964255 0 0 + 0.0952885 0 0.03 + 0.0935737 0 0.03 + 0.091859 0 0.03 + 0.0901442 0 0.03 + 0.0884295 0 0.03 + 0.0867147 0 0.03 + 0.1 0 -0.00357452 + 0.1 0 -0.00428206 + 0.1 0 -0.00515926 + 0.1 0 -0.00628309 + 0.1 0 -0.00776959 + 0.1 0 -0.00953435 + 0.1 0 -0.0113562 + 0.1 0 -0.0131781 + 0.102727 0 -0.015 + 0.105455 0 -0.015 + 0.108182 0 -0.015 + 0.110909 0 -0.015 + 0.113636 0 -0.015 + 0.116364 0 -0.015 + 0.119091 0 -0.015 + 0.121818 0 -0.015 + 0.124545 0 -0.015 + 0.127273 0 -0.015 + 0.13 0 -0.0132853 + 0.13 0 -0.0115705 + 0.13 0 -0.00985577 + 0.13 0 -0.00814103 + 0.13 0 -0.00642629 + 0.13 0 -0.00471154 + 0.085 0 0.03 + 0.085 0 0 + 0.1 0 -0.015 + 0.13 0 -0.015 + + + 598 300 301 + 599 329 330 + 326 327 600 + 335 336 601 + 604 603 602 + 367 368 605 + 370 371 606 + 899 779 778 + 785 784 908 + 901 781 780 + 783 782 906 + 788 787 786 + 147 789 148 + 893 655 206 + 198 199 894 + 898 897 896 + 931 886 885 + 892 891 940 + 933 888 887 + 890 889 938 + 197 198 895 + + + 0 21 1 22 + 1 22 2 23 + 2 23 3 24 + 3 24 4 25 + 4 25 5 26 + 5 26 6 27 + 6 27 7 28 + 7 28 8 29 + 8 29 9 30 + 9 30 10 31 + 10 31 11 32 + 11 32 12 33 + 12 33 13 34 + 13 34 14 35 + 14 35 15 36 + 15 36 16 37 + 16 37 17 38 + 17 38 18 39 + 18 39 19 40 + 19 40 20 41 + 21 42 22 43 + 22 43 23 44 + 23 44 24 45 + 24 45 25 46 + 25 46 26 47 + 26 47 27 48 + 27 48 28 49 + 28 49 29 50 + 29 50 30 51 + 30 51 31 52 + 31 52 32 53 + 32 53 33 54 + 33 54 34 55 + 34 55 35 56 + 35 56 36 57 + 36 57 37 58 + 37 58 38 59 + 38 59 39 60 + 39 60 40 61 + 40 61 41 62 + 42 63 43 64 + 43 64 44 65 + 44 65 45 66 + 45 66 46 67 + 46 67 47 68 + 47 68 48 69 + 48 69 49 70 + 49 70 50 71 + 50 71 51 72 + 51 72 52 73 + 52 73 53 74 + 53 74 54 75 + 54 75 55 76 + 55 76 56 77 + 56 77 57 78 + 57 78 58 79 + 58 79 59 80 + 59 80 60 81 + 60 81 61 82 + 61 82 62 83 + 63 84 64 85 + 64 85 65 86 + 65 86 66 87 + 66 87 67 88 + 67 88 68 89 + 68 89 69 90 + 69 90 70 91 + 70 91 71 92 + 71 92 72 93 + 72 93 73 94 + 73 94 74 95 + 74 95 75 96 + 75 96 76 97 + 76 97 77 98 + 77 98 78 99 + 78 99 79 100 + 79 100 80 101 + 80 101 81 102 + 81 102 82 103 + 82 103 83 104 + 84 105 85 106 + 85 106 86 107 + 86 107 87 108 + 87 108 88 109 + 88 109 89 110 + 89 110 90 111 + 90 111 91 112 + 91 112 92 113 + 92 113 93 114 + 93 114 94 115 + 94 115 95 116 + 95 116 96 117 + 96 117 97 118 + 97 118 98 119 + 98 119 99 120 + 99 120 100 121 + 100 121 101 122 + 101 122 102 123 + 102 123 103 124 + 103 124 104 125 + 105 126 106 127 + 106 127 107 128 + 107 128 108 129 + 108 129 109 130 + 109 130 110 131 + 110 131 111 132 + 111 132 112 133 + 112 133 113 134 + 113 134 114 135 + 114 135 115 136 + 115 136 116 137 + 116 137 117 138 + 117 138 118 139 + 118 139 119 140 + 119 140 120 141 + 120 141 121 142 + 121 142 122 143 + 122 143 123 144 + 123 144 124 145 + 124 145 125 146 + 205 202 174 0 + 202 201 0 21 + 201 200 21 42 + 200 199 42 63 + 199 198 63 84 + 198 197 84 105 + 197 196 105 126 + 196 206 126 175 + 126 175 127 176 + 127 176 128 177 + 128 177 129 178 + 129 178 130 179 + 130 179 131 180 + 131 180 132 181 + 132 181 133 182 + 133 182 134 183 + 134 183 135 184 + 135 184 136 185 + 136 185 137 186 + 137 186 138 187 + 138 187 139 188 + 139 188 140 189 + 140 189 141 190 + 141 190 142 191 + 142 191 143 192 + 143 192 144 193 + 144 193 145 194 + 145 194 146 195 + 146 195 147 204 + 125 146 148 147 + 104 125 149 148 + 83 104 150 149 + 62 83 151 150 + 41 62 152 151 + 20 41 153 152 + 154 20 203 153 + 155 19 154 20 + 156 18 155 19 + 157 17 156 18 + 158 16 157 17 + 159 15 158 16 + 160 14 159 15 + 161 13 160 14 + 162 12 161 13 + 163 11 162 12 + 164 10 163 11 + 165 9 164 10 + 166 8 165 9 + 167 7 166 8 + 168 6 167 7 + 169 5 168 6 + 170 4 169 5 + 171 3 170 4 + 172 2 171 3 + 173 1 172 2 + 174 0 173 1 + 207 678 208 623 + 208 623 209 624 + 209 624 210 625 + 210 625 211 626 + 211 626 212 627 + 212 627 213 628 + 213 628 214 629 + 214 629 215 630 + 215 630 216 631 + 216 631 217 632 + 217 632 218 633 + 218 633 219 634 + 219 634 220 635 + 220 635 221 636 + 221 636 222 637 + 222 637 223 638 + 223 638 224 639 + 224 639 607 677 + 224 607 225 608 + 225 608 226 609 + 226 609 227 610 + 227 610 228 611 + 228 611 229 612 + 229 612 230 613 + 230 613 231 614 + 231 614 232 615 + 232 615 233 616 + 233 616 234 617 + 234 617 235 618 + 235 618 236 619 + 236 619 237 620 + 237 620 238 621 + 238 621 239 622 + 239 622 195 204 + 239 195 240 194 + 240 194 241 193 + 241 193 242 192 + 242 192 243 191 + 243 191 244 190 + 244 190 245 189 + 245 189 246 188 + 246 188 247 187 + 247 187 248 186 + 248 186 249 185 + 249 185 250 184 + 250 184 251 183 + 251 183 252 182 + 252 182 253 181 + 253 181 254 180 + 254 180 255 179 + 255 179 256 178 + 256 178 257 177 + 257 177 258 176 + 258 176 259 175 + 259 175 655 206 + 259 655 260 654 + 260 654 261 653 + 261 653 262 652 + 262 652 263 651 + 263 651 264 650 + 264 650 265 649 + 265 649 266 648 + 266 648 267 647 + 267 647 268 646 + 268 646 269 645 + 269 645 270 644 + 270 644 271 643 + 271 643 272 642 + 272 642 273 641 + 273 641 274 640 + 274 640 660 679 + 274 660 275 661 + 275 661 276 662 + 276 662 277 663 + 277 663 278 664 + 278 664 279 665 + 279 665 280 666 + 280 666 281 667 + 281 667 282 668 + 282 668 283 669 + 283 669 284 670 + 284 670 285 671 + 285 671 286 672 + 286 672 287 673 + 287 673 288 674 + 288 674 289 675 + 289 675 290 676 + 290 676 291 680 + 291 680 292 656 + 292 656 293 657 + 293 657 294 658 + 294 658 295 659 + 295 659 207 678 + 296 207 297 208 + 297 208 298 209 + 298 209 299 210 + 299 210 300 211 + 300 211 301 212 + 301 212 302 213 + 302 213 303 214 + 303 214 304 215 + 304 215 305 216 + 305 216 306 217 + 306 217 307 218 + 307 218 308 219 + 308 219 309 220 + 309 220 310 221 + 310 221 311 222 + 311 222 312 223 + 312 223 225 224 + 312 225 313 226 + 313 226 314 227 + 314 227 315 228 + 315 228 316 229 + 316 229 317 230 + 317 230 318 231 + 318 231 319 232 + 319 232 320 233 + 320 233 321 234 + 321 234 322 235 + 322 235 323 236 + 323 236 324 237 + 324 237 325 238 + 325 238 240 239 + 325 240 326 241 + 326 241 327 242 + 327 242 328 243 + 328 243 329 244 + 329 244 330 245 + 330 245 331 246 + 331 246 332 247 + 332 247 333 248 + 333 248 334 249 + 334 249 335 250 + 335 250 336 251 + 336 251 337 252 + 337 252 338 253 + 338 253 339 254 + 339 254 340 255 + 340 255 341 256 + 341 256 342 257 + 342 257 343 258 + 343 258 260 259 + 343 260 344 261 + 344 261 345 262 + 345 262 346 263 + 346 263 347 264 + 347 264 348 265 + 348 265 349 266 + 349 266 350 267 + 350 267 351 268 + 351 268 352 269 + 352 269 353 270 + 353 270 354 271 + 354 271 355 272 + 355 272 356 273 + 356 273 275 274 + 356 275 357 276 + 357 276 358 277 + 358 277 359 278 + 359 278 360 279 + 360 279 361 280 + 361 280 362 281 + 362 281 363 282 + 363 282 364 283 + 364 283 365 284 + 365 284 366 285 + 366 285 367 286 + 367 286 368 287 + 368 287 369 288 + 369 288 370 289 + 370 289 371 290 + 371 290 372 291 + 372 291 373 292 + 373 292 374 293 + 374 293 375 294 + 375 294 376 295 + 376 295 296 207 + 377 392 378 393 + 380 394 381 395 + 381 395 382 396 + 382 396 383 397 + 383 397 384 398 + 384 398 385 399 + 385 399 386 400 + 386 400 387 401 + 387 401 388 402 + 388 402 389 403 + 389 403 390 404 + 390 404 391 405 + 391 405 392 406 + 392 406 393 407 + 394 408 395 409 + 395 409 396 410 + 396 410 397 411 + 397 411 398 412 + 398 412 399 413 + 399 413 400 414 + 400 414 401 415 + 401 415 402 416 + 402 416 403 417 + 403 417 404 418 + 404 418 405 419 + 405 419 406 420 + 406 420 407 421 + 408 422 409 423 + 409 423 410 424 + 410 424 411 425 + 411 425 412 426 + 412 426 413 427 + 413 427 414 428 + 414 428 415 429 + 415 429 416 430 + 416 430 417 431 + 417 431 418 432 + 418 432 419 433 + 419 433 420 434 + 420 434 421 435 + 422 436 423 437 + 423 437 424 438 + 424 438 425 439 + 425 439 426 440 + 426 440 427 441 + 427 441 428 442 + 428 442 429 443 + 429 443 430 444 + 430 444 431 445 + 431 445 432 446 + 432 446 433 447 + 433 447 434 448 + 434 448 435 449 + 436 450 437 451 + 437 451 438 452 + 438 452 439 453 + 439 453 440 454 + 440 454 441 455 + 441 455 442 456 + 442 456 443 457 + 443 457 444 458 + 444 458 445 459 + 445 459 446 460 + 446 460 447 461 + 447 461 448 462 + 448 462 449 463 + 450 464 451 465 + 451 465 452 466 + 452 466 453 467 + 453 467 454 468 + 454 468 455 469 + 455 469 456 470 + 456 470 457 471 + 457 471 458 472 + 458 472 459 473 + 459 473 460 474 + 460 474 461 475 + 461 475 462 476 + 462 476 463 477 + 464 478 465 479 + 465 479 466 480 + 466 480 467 481 + 467 481 468 482 + 468 482 469 483 + 469 483 470 484 + 470 484 471 485 + 471 485 472 486 + 472 486 473 487 + 473 487 474 488 + 474 488 475 489 + 475 489 476 490 + 476 490 477 491 + 478 492 479 493 + 479 493 480 494 + 480 494 481 495 + 481 495 482 496 + 482 496 483 497 + 483 497 484 498 + 484 498 485 499 + 485 499 486 500 + 486 500 487 501 + 487 501 488 502 + 488 502 489 503 + 489 503 490 504 + 490 504 491 505 + 492 506 493 507 + 493 507 494 508 + 494 508 495 509 + 495 509 496 510 + 496 510 497 511 + 497 511 498 512 + 498 512 499 513 + 499 513 500 514 + 500 514 501 515 + 501 515 502 516 + 502 516 503 517 + 503 517 504 518 + 504 518 505 519 + 506 520 507 521 + 507 521 508 522 + 508 522 509 523 + 509 523 510 524 + 510 524 511 525 + 511 525 512 526 + 512 526 513 527 + 513 527 514 528 + 514 528 515 529 + 515 529 516 530 + 516 530 517 531 + 517 531 518 532 + 518 532 519 533 + 520 534 521 535 + 521 535 522 536 + 522 536 523 537 + 523 537 524 538 + 524 538 525 539 + 525 539 526 540 + 526 540 527 541 + 527 541 528 542 + 528 542 529 543 + 529 543 530 544 + 530 544 531 545 + 531 545 532 546 + 532 546 533 547 + 535 548 536 549 + 536 549 537 550 + 537 550 538 551 + 538 551 539 552 + 539 552 540 553 + 540 553 541 554 + 541 554 542 555 + 542 555 543 556 + 543 556 544 557 + 544 557 545 558 + 545 558 546 559 + 546 559 547 560 + 549 561 550 562 + 550 562 551 563 + 551 563 552 564 + 552 564 553 565 + 553 565 554 566 + 554 566 555 567 + 555 567 556 568 + 556 568 557 569 + 557 569 558 570 + 558 570 559 571 + 559 571 560 572 + 384 385 307 306 + 307 308 384 383 + 308 309 383 382 + 309 310 382 381 + 310 311 381 380 + 311 312 380 313 + 303 304 388 387 + 304 305 387 386 + 305 306 386 385 + 301 302 390 389 + 302 303 389 388 + 598 301 391 390 + 392 377 391 598 + 377 299 598 300 + 574 576 520 534 + 576 599 534 573 + 599 575 329 328 + 599 330 573 331 + 464 450 319 318 + 450 436 318 317 + 436 422 317 316 + 422 408 316 315 + 408 394 315 314 + 394 380 314 313 + 327 328 600 575 + 326 600 325 324 + 324 600 323 578 + 600 575 578 574 + 574 520 578 579 + 520 506 579 577 + 478 464 320 319 + 492 478 321 320 + 321 577 492 506 + 333 334 580 548 + 581 549 601 548 + 581 601 337 336 + 604 581 338 337 + 335 601 334 548 + 349 350 566 567 + 350 351 567 568 + 351 352 568 569 + 352 353 569 570 + 353 354 570 571 + 354 355 571 572 + 603 583 602 584 + 563 562 582 585 + 348 349 565 566 + 347 348 564 565 + 339 603 338 604 + 370 606 369 586 + 606 393 586 407 + 367 605 435 421 + 369 586 368 605 + 586 407 605 421 + 435 449 367 366 + 449 463 366 365 + 463 477 365 364 + 477 491 364 363 + 491 505 363 362 + 505 519 362 361 + 519 533 361 360 + 533 547 360 359 + 547 560 359 358 + 560 572 358 357 + 572 355 357 356 + 374 375 379 376 + 298 299 378 377 + 575 599 574 576 + 579 577 322 321 + 579 322 578 323 + 580 548 587 535 + 587 535 573 534 + 333 580 332 587 + 332 587 331 573 + 581 604 549 561 + 583 603 340 339 + 345 582 584 585 + 561 604 588 602 + 602 584 588 585 + 585 562 588 561 + 341 589 340 583 + 341 342 589 590 + 342 343 590 344 + 590 344 584 345 + 589 590 583 584 + 296 591 376 379 + 296 297 591 592 + 297 298 592 378 + 346 347 593 564 + 564 563 593 582 + 582 345 593 346 + 606 371 594 372 + 594 372 595 373 + 595 373 379 374 + 592 378 596 393 + 393 606 596 594 + 596 594 597 595 + 597 595 591 379 + 591 592 597 596 + 778 779 681 686 + 681 686 682 687 + 682 687 683 688 + 683 688 684 689 + 684 689 685 690 + 779 691 686 692 + 686 692 687 693 + 687 693 688 694 + 688 694 689 695 + 689 695 690 696 + 691 780 692 697 + 692 697 693 698 + 693 698 694 699 + 694 699 695 700 + 695 700 696 701 + 780 781 697 702 + 697 702 698 703 + 698 703 699 704 + 699 704 700 705 + 700 705 701 706 + 781 707 702 708 + 702 708 703 709 + 703 709 704 710 + 704 710 705 711 + 705 711 706 712 + 707 713 708 714 + 708 714 709 715 + 709 715 710 716 + 710 716 711 717 + 711 717 712 718 + 713 719 714 720 + 714 720 715 721 + 715 721 716 722 + 716 722 717 723 + 717 723 718 724 + 719 725 720 726 + 720 726 721 727 + 721 727 722 728 + 722 728 723 729 + 723 729 724 730 + 725 731 726 732 + 726 732 727 733 + 727 733 728 734 + 728 734 729 735 + 729 735 730 736 + 731 737 732 738 + 732 738 733 739 + 733 739 734 740 + 734 740 735 741 + 735 741 736 742 + 737 782 738 743 + 738 743 739 744 + 739 744 740 745 + 740 745 741 746 + 741 746 742 747 + 782 783 743 748 + 743 748 744 749 + 744 749 745 750 + 745 750 746 751 + 746 751 747 752 + 783 753 748 754 + 748 754 749 755 + 749 755 750 756 + 750 756 751 757 + 751 757 752 758 + 753 784 754 759 + 754 759 755 760 + 755 760 756 761 + 756 761 757 762 + 757 762 758 763 + 784 785 759 764 + 759 764 760 765 + 760 765 761 766 + 761 766 762 786 + 762 786 763 787 + 768 769 767 787 + 607 677 685 917 + 608 607 690 685 + 696 701 609 610 + 609 608 696 690 + 763 771 758 770 + 622 621 771 770 + 772 770 620 621 + 752 772 619 620 + 611 610 706 701 + 706 712 611 612 + 613 612 718 712 + 718 724 613 614 + 747 752 618 619 + 615 614 730 724 + 616 615 736 730 + 742 747 617 618 + 617 616 742 736 + 684 685 918 917 + 918 919 684 683 + 682 683 920 919 + 921 922 681 778 + 681 682 921 920 + 899 900 779 691 + 899 778 947 922 + 784 753 908 907 + 785 908 909 948 + 901 902 781 707 + 901 780 900 691 + 782 737 906 905 + 783 906 753 907 + 903 713 902 707 + 731 904 737 905 + 764 785 910 909 + 765 764 911 910 + 912 913 766 786 + 766 765 912 911 + 153 773 203 916 + 788 786 914 913 + 774 767 788 787 + 915 916 774 773 + 774 788 915 914 + 789 763 769 787 + 149 769 150 768 + 771 763 775 789 + 775 789 204 147 + 775 204 771 622 + 772 752 770 758 + 903 904 725 731 + 903 725 713 719 + 773 153 776 152 + 776 767 773 774 + 789 769 148 149 + 776 152 777 151 + 151 150 777 768 + 768 767 777 776 + 897 898 790 805 + 790 805 791 806 + 791 806 792 807 + 792 807 793 808 + 793 808 794 809 + 794 809 795 810 + 795 810 796 811 + 796 811 797 812 + 797 812 798 813 + 798 813 799 814 + 799 814 800 815 + 800 815 801 816 + 801 816 802 817 + 802 817 803 818 + 803 818 804 819 + 898 820 805 821 + 805 821 806 822 + 806 822 807 823 + 807 823 808 824 + 808 824 809 825 + 809 825 810 826 + 810 826 811 827 + 811 827 812 828 + 812 828 813 829 + 813 829 814 830 + 814 830 815 831 + 815 831 816 832 + 816 832 817 833 + 817 833 818 834 + 818 834 819 835 + 820 836 821 837 + 821 837 822 838 + 822 838 823 839 + 823 839 824 840 + 824 840 825 841 + 825 841 826 842 + 826 842 827 843 + 827 843 828 844 + 828 844 829 845 + 829 845 830 846 + 830 846 831 847 + 831 847 832 848 + 832 848 833 849 + 833 849 834 850 + 834 850 835 851 + 836 852 837 853 + 837 853 838 854 + 838 854 839 855 + 839 855 840 856 + 840 856 841 857 + 841 857 842 858 + 842 858 843 859 + 843 859 844 860 + 844 860 845 861 + 845 861 846 862 + 846 862 847 863 + 847 863 848 864 + 848 864 849 865 + 849 865 850 866 + 850 866 851 867 + 852 885 853 886 + 853 886 854 868 + 854 868 855 887 + 855 887 856 888 + 856 888 857 869 + 857 869 858 870 + 858 870 859 871 + 859 871 860 872 + 860 872 861 873 + 861 873 862 874 + 862 874 863 889 + 863 889 864 890 + 864 890 865 875 + 865 875 866 891 + 866 891 867 892 + 803 804 641 640 + 642 643 802 801 + 802 803 642 641 + 893 790 877 791 + 893 877 655 654 + 877 878 654 653 + 878 792 653 652 + 800 801 644 643 + 644 645 800 799 + 798 799 646 645 + 646 647 798 797 + 651 652 793 792 + 796 797 648 647 + 795 796 649 648 + 650 651 794 793 + 794 795 650 649 + 880 894 897 879 + 199 200 894 879 + 929 930 852 885 + 928 929 836 852 + 820 898 927 926 + 927 928 820 836 + 881 879 201 200 + 898 896 926 925 + 879 882 897 896 + 924 925 882 896 + 931 932 886 868 + 931 885 949 930 + 891 875 940 939 + 892 940 941 950 + 933 934 888 869 + 933 887 932 868 + 889 874 938 937 + 890 938 875 939 + 935 870 934 869 + 873 936 874 937 + 819 835 945 944 + 946 679 804 640 + 804 819 946 945 + 943 944 851 835 + 867 892 942 941 + 942 943 867 851 + 198 894 895 880 + 880 897 876 790 + 876 790 883 893 + 883 893 196 206 + 880 876 895 883 + 895 883 197 196 + 791 792 877 878 + 881 201 884 202 + 884 202 923 205 + 879 881 882 884 + 882 884 924 923 + 935 936 872 873 + 935 872 870 871 + + diff --git a/src/SIM/SIM1D.C b/src/SIM/SIM1D.C index 203258204..ce699a0de 100644 --- a/src/SIM/SIM1D.C +++ b/src/SIM/SIM1D.C @@ -148,7 +148,7 @@ bool SIM1D::parseGeometryTag (const tinyxml2::XMLElement* elem) utl::getAttribute(elem,"offset",offset); const tinyxml2::XMLElement* child = elem->FirstChildElement("connection"); - for (; child; child = child->NextSiblingElement()) + for (; child; child = child->NextSiblingElement("connection")) { ASM::Interface ifc; if (utl::getAttribute(child,"master",ifc.master)) diff --git a/src/SIM/SIM2D.C b/src/SIM/SIM2D.C index d1094e48e..51d09171c 100644 --- a/src/SIM/SIM2D.C +++ b/src/SIM/SIM2D.C @@ -187,7 +187,7 @@ bool SIM2D::parseGeometryTag (const tinyxml2::XMLElement* elem) std::vector top; const tinyxml2::XMLElement* child = elem->FirstChildElement("connection"); - for (; child; child = child->NextSiblingElement()) + for (; child; child = child->NextSiblingElement("connection")) { ASM::Interface ifc; bool rever = false, periodic = false; diff --git a/src/SIM/SIM3D.C b/src/SIM/SIM3D.C index e37a6f1b5..fdce6c69c 100644 --- a/src/SIM/SIM3D.C +++ b/src/SIM/SIM3D.C @@ -178,7 +178,7 @@ bool SIM3D::parseGeometryTag (const tinyxml2::XMLElement* elem) utl::getAttribute(elem,"offset",offset); const tinyxml2::XMLElement* child = elem->FirstChildElement("connection"); - for (; child; child = child->NextSiblingElement()) + for (; child; child = child->NextSiblingElement("connection")) { ASM::Interface ifc; bool periodic = false; diff --git a/src/SIM/SIMinput.C b/src/SIM/SIMinput.C index 4fe9da008..d1e1727d3 100644 --- a/src/SIM/SIMinput.C +++ b/src/SIM/SIMinput.C @@ -451,7 +451,7 @@ bool SIMinput::parseBCTag (const tinyxml2::XMLElement* elem) else if (!strcasecmp(elem->Value(),"propertycodes")) { const tinyxml2::XMLElement* code = elem->FirstChildElement("code"); - for (; code; code = code->NextSiblingElement()) + for (; code; code = code->NextSiblingElement("code")) { int icode = 0; utl::getAttribute(code,"value",icode); diff --git a/src/SIM/SIMoutput.C b/src/SIM/SIMoutput.C index a5988d162..e9eb1ae57 100644 --- a/src/SIM/SIMoutput.C +++ b/src/SIM/SIMoutput.C @@ -106,8 +106,7 @@ bool SIMoutput::parseOutputTag (const tinyxml2::XMLElement* elem) { IFEM::cout <<" Parsing <"<< elem->Value() <<">"<< std::endl; - const char* funcval = utl::getValue(elem,"function"); - if (funcval) + if (const char* funcval = utl::getValue(elem,"function"); funcval) { std::string name; utl::getAttribute(elem,"name",name); @@ -138,10 +137,13 @@ bool SIMoutput::parseOutputTag (const tinyxml2::XMLElement* elem) points.back().second.emplace_back(newPoint); }; + Vec3 X0; // Optional coordinate shift for Cartesian coordinates + utl::getAttribute(elem,"shift",X0); + // Parse the result point specifications. // Can either be explicit points, lines or a grid. const tinyxml2::XMLElement* point = elem->FirstChildElement("point"); - for (int i = 1; point; i++, point = point->NextSiblingElement()) + for (int i = 1; point; i++, point = point->NextSiblingElement("point")) { ResultPoint thePoint(3); if (int patch = 0; utl::getAttribute(point,"patch",patch) && patch > 0) @@ -150,21 +152,31 @@ bool SIMoutput::parseOutputTag (const tinyxml2::XMLElement* elem) IFEM::cout <<"\tPoint "<< i <<": P"<< thePoint.patch; if (utl::getAttribute(point,"node",thePoint.inod) && thePoint.inod > 0) IFEM::cout <<" node = "<< thePoint.inod; + else if (point->FirstChild() && + utl::parseVec(thePoint.X,point->FirstChild()->Value())) + { + // Cartesian coordinates of a spatial point is specified, + thePoint.iel = -1; // need to search for the matching element/node + thePoint.X += X0; + IFEM::cout <<" X = "<< thePoint.X; + } else + { IFEM::cout <<" xi ="; - if (utl::getAttribute(point,"u",thePoint.u[0])) - IFEM::cout <<' '<< thePoint.u[0]; - if (utl::getAttribute(point,"v",thePoint.u[1])) - IFEM::cout <<' '<< thePoint.u[1]; - if (utl::getAttribute(point,"w",thePoint.u[2])) - IFEM::cout <<' '<< thePoint.u[2]; + if (utl::getAttribute(point,"u",thePoint.u[0])) + IFEM::cout <<' '<< thePoint.u[0]; + if (utl::getAttribute(point,"v",thePoint.u[1])) + IFEM::cout <<' '<< thePoint.u[1]; + if (utl::getAttribute(point,"w",thePoint.u[2])) + IFEM::cout <<' '<< thePoint.u[2]; + } IFEM::cout << std::endl; addPoint(thePoint); } const tinyxml2::XMLElement* line = elem->FirstChildElement("line"); - for (int j = 1; line; j++, line = line->NextSiblingElement()) + for (int j = 1; line; j++, line = line->NextSiblingElement("line")) { ResultPoint thePoint(3); if (int patch = 0; utl::getAttribute(line,"patch",patch) && patch > 0) @@ -206,26 +218,22 @@ bool SIMoutput::parseOutputTag (const tinyxml2::XMLElement* elem) } for (const tinyxml2::XMLElement* child = elem->FirstChildElement("elements"); - child; child = child->NextSiblingElement()) + child; child = child->NextSiblingElement("elements")) if (const tinyxml2::XMLNode* elist = child->FirstChild(); elist) { ResultPoint thePoint; if (int patch = 0; utl::getAttribute(child,"patch",patch) && patch > 0) thePoint.patch = patch; - IntVec elements; - utl::parseIntegers(elements,elist->Value()); - for (int iel : elements) - { - thePoint.iel = iel; - addPoint(thePoint); - } - - if (!elements.empty()) + if (IntVec elements; utl::parseIntegers(elements,elist->Value())) { IFEM::cout <<"\tElement center output: P"<< thePoint.patch <<" ("; - for (const ResultPoint& rp : myPoints.back().second) - IFEM::cout <<" "<< rp.iel; + for (int iel : elements) + { + IFEM::cout <<" "<< iel; + thePoint.iel = iel; + addPoint(thePoint); + } IFEM::cout <<" )"<< std::endl; } } @@ -236,7 +244,7 @@ bool SIMoutput::parseOutputTag (const tinyxml2::XMLElement* elem) else if (grid) idxGrid = myPoints.size(); - for (int g = 1; grid; g++, grid = grid->NextSiblingElement()) + for (int g = 1; grid; g++, grid = grid->NextSiblingElement("grid")) { ResultPoint thePoint; if (int patch = 0; utl::getAttribute(grid,"patch",patch) && patch > 0) @@ -413,7 +421,7 @@ void SIMoutput::preprocessResPtGroup (std::string& ptFile, ResPointVec& points) { if (points.empty()) return; - // Check if we have Cartesian grid input + // Check if we have a regular Cartesian grid input size_t nX = 0; IntMat pointMap; if (points.front().inod < 0 && points.back().inod < 0) @@ -426,23 +434,52 @@ void SIMoutput::preprocessResPtGroup (std::string& ptFile, ResPointVec& points) nX = 0; } + // Check if we have general spatial points input (Cartesian coordinates) + using PatchPoints = std::pair,Vec3Vec>; + std::map pchPoints; + if (pointMap.empty()) + for (ResultPoint& pt : points) + if (pt.iel < 0 || pt.inod < 0) + if (ASMbase* pch = this->getPatch(pt.patch,true); pch && !pch->empty()) + { + pchPoints[pch].second.emplace_back(pt.X); + pt.iel = -pchPoints[pch].second.size(); + } + + // Find element and/or local parameters of the spatial points + for (std::pair& pch : pchPoints) + if (!pch.first->findPoints(pch.second.second,pch.second.first)) + IFEM::cout <<" ** Failed to locate some result points."<< std::endl; + int iPoint = 0; size_t iX = 0, iY = 0; for (ResPointVec::iterator pit = points.begin(); pit != points.end();) { bool pointIsOK = false; - ASMbase* pch = this->getPatch(pit->patch,true); - if (pch && !pch->empty()) + if (ASMbase* pch = this->getPatch(pit->patch,true); pch && !pch->empty()) { if ((pointIsOK = pit->inod > 0)) // A nodal number is specified pit->X = pch->getCoord(pit->inod); - else if (pit->inod < 0) // A spatial point is specified - pointIsOK = fabs(pch->findPoint(pit->X,pit->u)) < 1.0e-3; - + else if (pit->iel < 0 && !pchPoints.empty()) + { + // A spatial point is specified, use the point searching results + size_t iP = -pit->iel-1; + std::map::const_iterator it = pchPoints.find(pch); + if (it != pchPoints.end() && iP < it->second.first.size()) + if ((pointIsOK = fabs(it->second.first[iP].dist) < 1.0e-3)) + { + pit->npar = pch->getNoParamDim(); + pit->X = it->second.second[iP]; + pit->iel = it->second.first[iP].iel; + memcpy(pit->u,it->second.first[iP].u,sizeof(double)*3); + } + } else if (pit->npar > 0) // A parametric point is specified + { + pit->npar = pch->getNoParamDim(); pointIsOK = (pit->inod = pch->evalPoint(pit->u,pit->u,pit->X)) >= 0; - + } else if (pit->iel < 1) // All nodes or element centers pointIsOK = true; @@ -455,8 +492,6 @@ void SIMoutput::preprocessResPtGroup (std::string& ptFile, ResPointVec& points) if (pointIsOK) { - if (pit->npar > 0) - pit->npar = pch->getNoParamDim(); ++pit; ++iPoint; if (nX > 0) @@ -498,10 +533,24 @@ void SIMoutput::preprocessResPtGroup (std::string& ptFile, ResPointVec& points) IFEM::cout <<", node #"<< pt.inod; else if (pt.iel > 0) IFEM::cout <<", element #"<< pt.iel; - if (pt.inod > 0 && myModel.size() > 1) + ASMbase* pch = this->getPatch(pt.patch,true); + if (pch) { - ASMbase* pch = this->getPatch(pt.patch,true); - IFEM::cout <<", global #"<< pch->getNodeID(pt.inod); + int globalId = 0; + if (pt.inod > 0) + { + if (int nodeId = pch->getNodeID(pt.inod); + nodeId != pt.inod || myModel.size() > 1) + globalId = nodeId; + } + else if (pt.iel > 0) + { + if (int elmId = pch->getElmID(pt.iel); + elmId != pt.iel || myModel.size() > 1) + globalId = elmId; + } + if (globalId > 0) + IFEM::cout <<", global #"<< globalId; } if (pt.npar > 0 || pt.iel > 0) IFEM::cout <<", X = "<< pt.X << std::endl; @@ -509,6 +558,8 @@ void SIMoutput::preprocessResPtGroup (std::string& ptFile, ResPointVec& points) IFEM::cout <<", element centers will be used."<< std::endl; else IFEM::cout <<", nodal points will be used."<< std::endl; + if (pch && pt.iel > 0) + pch->printElmInfo(pt.iel,myProblem); } } @@ -539,7 +590,7 @@ void SIMoutput::preprocessResPtGroup (std::string& ptFile, ResPointVec& points) for (const ResultPoint& pt : points) { os << 0.0 <<" "; // dummy time - for (short int k = 0; k < pt.npar; k++) + for (short int k = 0; k < nsd; k++) os << std::setw(flWidth) << pt.X[k]; os << std::endl; } @@ -561,8 +612,7 @@ bool SIMoutput::merge (SIMbase* that, const std::map* old2new, for (ResultPoint& rp : rptp.second) { if (++iPoint == 1) IFEM::cout <<'\n'; - ASMbase* pch = that->getPatch(rp.patch,true); - if (pch && rp.inod > 0) + if (ASMbase* pch = that->getPatch(rp.patch,true); pch && rp.inod > 0) IFEM::cout <<"Result point #"<< iPoint <<": patch #"<< rp.patch <<", local node "<< rp.inod <<" global node "<< pch->getNodeID(rp.inod) << std::endl; @@ -1796,11 +1846,11 @@ bool SIMoutput::dumpMatlabGrid (std::ostream& os, const std::string& name, os <<"]="<< name << std::endl; // Write out nodes for the specified topology sets - TopologySet::const_iterator tit; for (const std::string& setname : sets) { std::set nodeSet; - if ((tit = myEntitys.find(setname)) != myEntitys.end()) + if (TopologySet::const_iterator tit = myEntitys.find(setname); + tit != myEntitys.end()) for (const TopItem& top : tit->second) if (ASMbase* pch = this->getPatch(top.patch); pch) if (top.idim+1 == pch->getNoParamDim()) @@ -2021,7 +2071,7 @@ bool SIMoutput::dumpResults (const Vector& psol, double time, } bool newLine = false; - // Convenience lambda returning proper newline character before identifier. + // Convenience lambda returning proper newline char before identifier. auto&& newLineChr = [&newLine]() { return newLine ? "\n\t\t" : ":\t"; }; if (formatted && !sol1.empty()) @@ -2124,6 +2174,8 @@ bool SIMoutput::evalResults (const Vectors& psol, const ResPointVec& gPoints, else if (pt.iel > 0) { elms.push_back(pt.iel); + for (short int k = 0; k < pt.npar; k++) + params[k].push_back(pt.u[k]); } else if (pt.npar < 0) { @@ -2226,7 +2278,8 @@ bool SIMoutput::evalResults (const Vectors& psol, const ResPointVec& gPoints, if (elms.empty()) ok = patch->evalSolution(tmp,*myProblem,params.data(),points.empty()); else // evaluate at specified elements - ok = patch->evalSolution(tmp,*myProblem,elms); + ok = patch->evalSolution(tmp,*myProblem,elms, + params.front().empty() ? nullptr : params.data()); if ((ok &= augment(sol2,tmp)) && getNames) for (size_t i = 0; i < myProblem->getNoFields(2); i++) diff --git a/src/SIM/TimeStep.C b/src/SIM/TimeStep.C index 4df2b38b9..96f69aeda 100644 --- a/src/SIM/TimeStep.C +++ b/src/SIM/TimeStep.C @@ -118,7 +118,7 @@ bool TimeStep::parse (const tinyxml2::XMLElement* elem) if (f2 > 1.0) f2 = 1.0; const tinyxml2::XMLElement* child = elem->FirstChildElement("step"); - for (; child; child = child->NextSiblingElement()) + for (; child; child = child->NextSiblingElement("step")) { double start = 0.0, end = 0.0; utl::getAttribute(child,"start",start); diff --git a/src/Utility/Test/TestUtilities.C b/src/Utility/Test/TestUtilities.C index d5fb5bd68..c28893b0a 100644 --- a/src/Utility/Test/TestUtilities.C +++ b/src/Utility/Test/TestUtilities.C @@ -10,38 +10,45 @@ //! //============================================================================== +#include "Catch2Support.h" #include "Utilities.h" #include "Vec3.h" #include "tinyxml2.h" -#include "Catch2Support.h" - #include -#include TEST_CASE("TestUtilities.ParseIntegers") { - const char* simple = "1"; - const char* ranged = "1:5"; - const char* multip = "1 3 5:8 10"; std::vector values1, values2, values3; - utl::parseIntegers(values1, simple); - utl::parseIntegers(values2, ranged); - utl::parseIntegers(values3, multip); + REQUIRE(utl::parseIntegers(values1,"1")); + REQUIRE(utl::parseIntegers(values2,"1:5")); + REQUIRE(utl::parseIntegers(values3,"1 3 5:8 10")); - int i; REQUIRE(values1.size() == 1); REQUIRE(values2.size() == 5); REQUIRE(values3.size() == 7); REQUIRE(values1.front() == 1); - for (i = 0; i < 5; i++) - REQUIRE(values2[i] == 1+i); - for (i = 0; i < 3; i++) + for (size_t i = 0; i < values2.size(); i++) + REQUIRE(values2[i] == 1+static_cast(i)); + for (int i = 0; i < 3; i++) REQUIRE(values3[i] == 1+2*i); - for (i = 2; i < 6; i++) + for (int i = 2; i < 6; i++) REQUIRE(values3[i] == 3+i); - REQUIRE(values3[6] == 10); + REQUIRE(values3.back() == 10); +} + + +TEST_CASE("TestUtilities.ParseVec3") +{ + Vec3 v1, v2, v3; + REQUIRE(utl::parseVec(v1,"1.2")); + REQUIRE(utl::parseVec(v2,"1 2")); + REQUIRE(utl::parseVec(v3,"1.2 3.4 5.6")); + + REQUIRE(v1.equal({1.2,0.0,0.0},0.0)); + REQUIRE(v2.equal({1.0,2.0,0.0},0.0)); + REQUIRE(v3.equal({1.2,3.4,5.6},0.0)); } @@ -50,7 +57,7 @@ TEST_CASE("TestUtilities.ParseKnots") std::string simple("xi 0 0.5 0.9 1.0"); strtok(const_cast(simple.c_str())," "); std::vector xi; - utl::parseKnots(xi); + REQUIRE(utl::parseKnots(xi)); REQUIRE(xi.size() == 4); REQUIRE_THAT(xi[0], WithinAbs(0.0, 1e-14)); REQUIRE_THAT(xi[1], WithinRel(0.5)); @@ -60,7 +67,7 @@ TEST_CASE("TestUtilities.ParseKnots") std::string graded5("xi C5 79 0.01 2.0"); strtok(const_cast(graded5.c_str())," "); xi.clear(); - utl::parseKnots(xi); + REQUIRE(utl::parseKnots(xi)); REQUIRE(xi.size() == 79); REQUIRE_THAT(xi[39], WithinRel(0.5)); REQUIRE_THAT(xi.front()+xi.back(), WithinRel(1.0)); @@ -118,29 +125,30 @@ TEST_CASE("TestUtilities.GetAttribute") tinyxml2::XMLDocument doc; doc.Parse(input); - REQUIRE(doc.RootElement() != nullptr); - const char* bTrue[4] = { "booltrue" , "boolone" , "boolon" , "boolyes" }; - const char* bFalse[4] = { "boolfalse", "boolzero", "booloff", "boolno" }; + const tinyxml2::XMLElement* root = doc.RootElement(); + REQUIRE(root != nullptr); - int i; - for (i = 0; i < 4; i++) { - const tinyxml2::XMLElement* elem = doc.RootElement()->FirstChildElement(bTrue[i]); + const tinyxml2::XMLElement* elem; + for (const char* bTrue : { "booltrue", "boolone", "boolon", "boolyes" }) + { + elem = root->FirstChildElement(bTrue); REQUIRE(elem != nullptr); bool b = false; REQUIRE(utl::getAttribute(elem, "bar", b)); REQUIRE(b); } - for (i = 0; i < 4; i++) { - const tinyxml2::XMLElement* elem = doc.RootElement()->FirstChildElement(bFalse[i]); + for (const char* bFalse : { "boolfalse", "boolzero", "booloff", "boolno" }) + { + elem = root->FirstChildElement(bFalse); REQUIRE(elem != nullptr); bool b = false; REQUIRE(utl::getAttribute(elem, "bar", b)); REQUIRE(!b); } - const tinyxml2::XMLElement* elem = doc.RootElement()->FirstChildElement("intval"); + elem = root->FirstChildElement("intval"); REQUIRE(elem != nullptr); int val1 = 0; REQUIRE(utl::getAttribute(elem, "bar", val1)); @@ -149,13 +157,13 @@ TEST_CASE("TestUtilities.GetAttribute") REQUIRE(utl::getAttribute(elem, "bar", val2)); REQUIRE(val2 == 1); - elem = doc.RootElement()->FirstChildElement("realval"); + elem = root->FirstChildElement("realval"); REQUIRE(elem != nullptr); double val3 = 0.0; REQUIRE(utl::getAttribute(elem, "bar", val3)); REQUIRE_THAT(val3, WithinRel(2.01)); - elem = doc.RootElement()->FirstChildElement("vecval1"); + elem = root->FirstChildElement("vecval1"); REQUIRE(elem != nullptr); Vec3 val4; REQUIRE(utl::getAttribute(elem, "bar", val4)); @@ -171,7 +179,7 @@ TEST_CASE("TestUtilities.GetAttribute") REQUIRE_THAT(val4.y, WithinRel(1.2)); REQUIRE_THAT(val4.z, WithinRel(1.2)); - elem = doc.RootElement()->FirstChildElement("vecval2"); + elem = root->FirstChildElement("vecval2"); REQUIRE(elem != nullptr); REQUIRE(utl::getAttribute(elem, "bar", val4)); REQUIRE_THAT(val4.x, WithinRel(1.2)); @@ -186,7 +194,7 @@ TEST_CASE("TestUtilities.GetAttribute") REQUIRE_THAT(val4.y, WithinRel(3.4)); REQUIRE_THAT(val4.z, WithinRel(3.4)); - elem = doc.RootElement()->FirstChildElement("vecval3"); + elem = root->FirstChildElement("vecval3"); REQUIRE(elem != nullptr); REQUIRE(utl::getAttribute(elem, "bar", val4)); REQUIRE_THAT(val4.x, WithinRel(1.2)); @@ -205,13 +213,8 @@ TEST_CASE("TestUtilities.GetAttribute") TEST_CASE("TestUtilities.GetDirs") { - int cmp1 = utl::getDirs(1); - int cmp2 = utl::getDirs(2); - int cmp3 = utl::getDirs(3); - int cmp4 = utl::getDirs(4); - - REQUIRE(cmp1 == 1); - REQUIRE(cmp2 == 12); - REQUIRE(cmp3 == 123); - REQUIRE(cmp4 == 1234); + REQUIRE(utl::getDirs(1) == 1); + REQUIRE(utl::getDirs(2) == 12); + REQUIRE(utl::getDirs(3) == 123); + REQUIRE(utl::getDirs(4) == 1234); } diff --git a/src/Utility/Utilities.C b/src/Utility/Utilities.C index 25baa530f..dac8ca72f 100644 --- a/src/Utility/Utilities.C +++ b/src/Utility/Utilities.C @@ -14,15 +14,16 @@ #include "Utilities.h" #include "Vec3.h" #include "tinyxml2.h" +#include #include #include -#include -void utl::parseIntegers (std::vector& values, const char* argv) +bool utl::parseIntegers (std::vector& values, const char* argv) { - if (!argv) return; + if (!argv) return false; + size_t old = values.size(); char* endp = const_cast(argv); while (endp && strlen(endp) > 0) { @@ -34,6 +35,23 @@ void utl::parseIntegers (std::vector& values, const char* argv) values.push_back(values.back()+1); } } + + return values.size() > old; +} + + +bool utl::parseVec (Vec3& value, const char* argv) +{ + if (!argv) return false; + + char* endp = const_cast(argv); + for (int i = 0; i < 3; i++) + if (endp && strlen(endp) > 0) + value[i] = strtod(endp,&endp); + else + value[i] = Real(0); + + return true; } diff --git a/src/Utility/Utilities.h b/src/Utility/Utilities.h index afda170b0..41b17745d 100644 --- a/src/Utility/Utilities.h +++ b/src/Utility/Utilities.h @@ -15,11 +15,9 @@ #define _UTILITIES_H #include "matrix.h" -#include -#include -#include #include #include +#include namespace tinyxml2 { class XMLElement; @@ -38,7 +36,12 @@ namespace utl //! \param[in] argv Character string with integer data //! //! \details An integer range is recognised through the syntax \a i:j. - void parseIntegers(std::vector& values, const char* argv); + bool parseIntegers(std::vector& values, const char* argv); + + //! \brief Parses a character string into a point vector. + //! \param[out] value The parsed point vector + //! \param[in] argv Character string with data + bool parseVec(Vec3& value, const char* argv); //! \brief Parses a (possibly graded) sequence of knot values. //! \param xi The knot value(s) is/are appended to this vector