From f3cdb440a2869e11accc048048650562253ceed9 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Mon, 15 Mar 2021 08:25:47 +0100 Subject: [PATCH] Fixed: Override getNoIntPoints and getNoBouPoints to reflect the respective integrate methods for LRSplines --- src/ASM/LR/ASMLRSpline.C | 89 +++++++++++++++++++++++++--------------- src/ASM/LR/ASMLRSpline.h | 13 +++--- src/ASM/LR/ASMu3D.C | 24 +++++++++++ src/ASM/LR/ASMu3D.h | 3 ++ 4 files changed, 92 insertions(+), 37 deletions(-) diff --git a/src/ASM/LR/ASMLRSpline.C b/src/ASM/LR/ASMLRSpline.C index e5baf8f5..94cc643d 100644 --- a/src/ASM/LR/ASMLRSpline.C +++ b/src/ASM/LR/ASMLRSpline.C @@ -99,24 +99,23 @@ void LR::generateThreadGroups (ThreadGroups& threadGroups, #ifdef USE_OPENMP if (omp_get_max_threads() > 1) { - for (int i = 0; i < 2; i++) - threadGroups[i].clear(); - IntMat& answer = threadGroups[0]; + threadGroups[0].clear(); + threadGroups[1].clear(); IntVec status(nElement,0); // status vector for elements: - std::vector> additionals; + std::vector additionals; if (!addConstraints.empty()) { additionals.resize(nElement); - for (auto e : lr->getAllElements()) { - for (const LRSpline* lr2 : addConstraints) { + for (LR::Element* e : lr->getAllElements()) { + for (LRSpline* lr2 : addConstraints) { int elB = lr2->getElementContaining(e->midpoint()); - for (auto b2 : lr2->getElement(elB)->support()) - for (auto el3 : b2->support()) { - std::vector midpoint = el3->midpoint(); - std::vector> points; + for (LR::Basisfunction* b2 : lr2->getElement(elB)->support()) + for (LR::Element* el3 : b2->support()) { + RealArray midpoint = el3->midpoint(); + std::vector points; if (lr2->nElements() != lr->nElements()) { - std::vector diff(midpoint.size()); + RealArray diff(midpoint.size()); for (size_t j = 0; j < midpoint.size(); ++j) diff[j] = (el3->getParmax(j) - el3->getParmin(j)) / 4.0; if (lr->dimension() == 2) @@ -136,7 +135,7 @@ void LR::generateThreadGroups (ThreadGroups& threadGroups, } else points = {midpoint}; - for (const std::vector& vec : points) + for (const RealArray& vec : points) additionals[e->getId()].insert(lr->getElementContaining(vec)); } } @@ -156,14 +155,14 @@ void LR::generateThreadGroups (ThreadGroups& threadGroups, // look for available elements IntVec thisColor; - for (auto e : lr->getAllElements()) { + for (LR::Element* e : lr->getAllElements()) { int i = e->getId(); if (status[i] == 0) { status[i] = nColors+1; thisColor.push_back(i); fixedElements++; - for (auto b : e->support()) // for all basisfunctions with support here - for (auto el2 : b->support()) {// for all elements this function supports + for (LR::Basisfunction* b : e->support()) + for (LR::Element* el2 : b->support()) { int j = el2->getId(); if (status[j] == 0) // if not assigned a color yet status[j] = -1; // set as unavailable (with current color) @@ -174,7 +173,7 @@ void LR::generateThreadGroups (ThreadGroups& threadGroups, } } } - answer.push_back(thisColor); + threadGroups[0].push_back(thisColor); } return; } @@ -472,30 +471,28 @@ Vec3 ASMLRSpline::getElementCenter (int iel) const bool ASMLRSpline::checkThreadGroups (const IntMat& groups, - const std::vector bases, + const std::vector& bases, const LR::LRSpline* threadBasis) { - size_t groupId = 1; bool ok = true; - for (const IntVec& elms : groups) { - size_t basisId = 1; - for (const LR::LRSpline* basis : bases) { - std::set nodes; - for (int elm : elms) { - std::vector midpoint = threadBasis->getElement(elm)->midpoint(); + for (size_t gId = 1; gId <= groups.size(); gId++) + for (size_t bId = 1; bId <= bases.size(); bId++) + { + IntSet nodes; + const LR::LRSpline* basis = bases[bId-1]; + for (int elm : groups[gId-1]) { + RealArray midpoint = threadBasis->getElement(elm)->midpoint(); int bElm = basis->getElementContaining(midpoint); - for (const LR::Basisfunction* func : basis->getElement(bElm)->support()) { - if (nodes.find(func->getId()) != nodes.end()) { - std::cerr << " ** Warning: Function " << func->getId() << " on basis " << basisId << " present for multiple elements in group " << groupId << std::endl; + for (const LR::Basisfunction* func : basis->getElement(bElm)->support()) + if (!nodes.insert(func->getId()).second) { + std::cerr <<" *** ASMLRSpline::checkThreadGroups: Function " + << func->getId() <<" on basis "<< bId + <<" is present for multiple elements in group "<< gId + << std::endl; ok = false; } - nodes.insert(func->getId()); - } } - ++basisId; } - ++groupId; - } return ok; } @@ -559,3 +556,31 @@ bool ASMLRSpline::getParameterDomain (Real2DMat& u, IntVec* corners) const return true; } + + +void ASMLRSpline::getNoIntPoints (size_t& nPt, size_t& nIPt) +{ + size_t nGp = 1; + if (nGauss > 0 && nGauss <= 10) + for (unsigned char d = 0; d < ndim; d++) + nGp *= nGauss; + else + { + // Use polynomial order to define number of quadrature points + int ng[3] = { 0, 0, 0 }; + int nG1 = 0; + this->getOrder(ng[0],ng[1],ng[2]); + for (unsigned char d = 0; d < ndim; d++) + nG1 = std::max(nG1,ng[d]+nGauss%10); + for (unsigned char d = 0; d < ndim; d++) + nGp = nG1*nG1; + } + + firstIp = nPt; + nPt += nel*nGp; + + // Count additional interface quadrature points + size_t nInterface = MLGE.size() - nel; + if (nInterface > 0 && nInterface != nel && nGauss > 0 && nGauss <= 10) + nIPt += nInterface*nGp/nGauss; +} diff --git a/src/ASM/LR/ASMLRSpline.h b/src/ASM/LR/ASMLRSpline.h index 8507c4b2..42d82472 100644 --- a/src/ASM/LR/ASMLRSpline.h +++ b/src/ASM/LR/ASMLRSpline.h @@ -178,6 +178,9 @@ public: //! \brief Returns the coordinates of the element center. virtual Vec3 getElementCenter(int iel) const; + //! \brief Computes the total number of integration points in this patch. + virtual void getNoIntPoints(size_t& nPt, size_t& nIPt); + protected: //! \brief Refines the mesh adaptively. //! \param[in] prm Input data used to control the mesh refinement @@ -188,15 +191,15 @@ protected: virtual int evalPoint(int iel, const double* param, Vec3& X) const = 0; //! \brief Santity check thread groups. - //! \param groups The generated thread groups - //! \param bases The bases to check for - //! \param threadBasis The LRSpline the element groups are derived from + //! \param[in] groups The generated thread groups + //! \param[in] bases The bases to check for + //! \param[in] threadBasis The LRSpline the element groups are derived from static bool checkThreadGroups(const IntMat& groups, - const std::vector bases, + const std::vector& bases, const LR::LRSpline* threadBasis); //! \brief Analyze and print thread group statistics. - //! \param groups The generated thread groups + //! \param[in] groups The generated thread groups static void analyzeThreadGroups(const IntMat& groups); LR::LRSpline* geo; //!< Pointer to the actual spline geometry object diff --git a/src/ASM/LR/ASMu3D.C b/src/ASM/LR/ASMu3D.C index 0b4d2e3c..8803a069 100644 --- a/src/ASM/LR/ASMu3D.C +++ b/src/ASM/LR/ASMu3D.C @@ -1982,6 +1982,30 @@ int ASMu3D::getCorner(int I, int J, int K, int basis) const } +void ASMu3D::getNoBouPoints (size_t& nPt, char ldim, char lindx) +{ + size_t nGp = 1; + if (nGauss > 0 && nGauss <= 10) + for (char d = 0; d < ldim; d++) + nGp *= nGauss; + else if (ldim == 2) + { + // Use polynomial order to define number of quadrature points + int p[3] = { 0, 0, 0 }; + this->getOrder(p[0],p[1],p[2]); + p[(lindx-1)/2] = 1; + int nG = std::max(std::max(p[0],p[1]),p[2]); + nGp = nG*nG; + } + else + nGp = 0; + + firstBp[lindx] = nPt; + + nPt += this->getNoBoundaryElms(lindx,ldim)*nGp; +} + + void ASMu3D::generateThreadGroups (const Integrand& integrand, bool silence, bool ignoreGlobalLM) { diff --git a/src/ASM/LR/ASMu3D.h b/src/ASM/LR/ASMu3D.h index 7b1cc60e..eef8a915 100644 --- a/src/ASM/LR/ASMu3D.h +++ b/src/ASM/LR/ASMu3D.h @@ -277,6 +277,9 @@ public: // These are the main computational methods of the ASM class hierarchy. // ==================================================================== + //! \brief Computes the number of boundary integration points in this patch. + virtual void getNoBouPoints(size_t& nPt, char ldim, char lindx); + //! \brief Evaluates an integral over the interior patch domain. //! \param integrand Object with problem-specific data and methods //! \param glbInt The integrated quantity