diff --git a/src/ASM/ASMunstruct.h b/src/ASM/ASMunstruct.h index d232f6bf7..47b3ae08f 100644 --- a/src/ASM/ASMunstruct.h +++ b/src/ASM/ASMunstruct.h @@ -46,6 +46,7 @@ namespace LR //! Utilities for LR-splines. IntVec options; //!< Parameters used to control the refinement IntVec elements; //!< 0-based indices of the elements to refine RealArray errors; //!< List of error indicators for the elements + std::vector MLGN; //!< MLGN mapping to use for multipatch //! \brief Default constructor. explicit RefineData(bool rs = false) : refShare(rs) {} diff --git a/src/ASM/LR/ASMLRSpline.C b/src/ASM/LR/ASMLRSpline.C index 67cc418a7..686bd609c 100644 --- a/src/ASM/LR/ASMLRSpline.C +++ b/src/ASM/LR/ASMLRSpline.C @@ -16,6 +16,7 @@ #include "LRSpline/LRSplineVolume.h" #include "ASMLRSpline.h" +#include "GlobalNodes.h" #include "Vec3.h" #include "Vec3Oper.h" #include "ThreadGroups.h" @@ -387,12 +388,12 @@ IntVec ASMLRSpline::getBoundaryCovered (const IntSet& nodes) const int numbEdges = (this->getNoParamDim() == 2) ? 4 : 6; for (int edge = 1; edge <= numbEdges; edge++) { - IntVec oneBoundary; // 1-based list of boundary nodes - this->getBoundaryNodes(edge,oneBoundary,1,1,0,true); - for (int i : nodes) - for (int j : oneBoundary) - if (geomB->getBasisfunction(i)->contains(*geomB->getBasisfunction(j-1))) - result.insert(j-1); + IntVec bnd = GlobalNodes::getBoundaryNodes(*refB, + this->getNoParamDim()-1, edge, 0); + for (const int i : nodes) + for (const int j : bnd) + if (refB->getBasisfunction(i)->contains(*refB->getBasisfunction(j))) + result.insert(j); } return IntVec(result.begin(), result.end()); @@ -404,7 +405,7 @@ IntVec ASMLRSpline::getOverlappingNodes (const IntSet& nodes, int dir) const IntSet result; for (int i : nodes) { - LR::Basisfunction* b = geomB->getBasisfunction(i); + const LR::Basisfunction* b = refB->getBasisfunction(i); for (LR::Element* el : b->support()) for (LR::Basisfunction* basis : el->support()) { diff --git a/src/ASM/LR/ASMu2D.C b/src/ASM/LR/ASMu2D.C index 649946642..88c3a53b1 100644 --- a/src/ASM/LR/ASMu2D.C +++ b/src/ASM/LR/ASMu2D.C @@ -21,6 +21,7 @@ #include "TimeDomain.h" #include "FiniteElement.h" #include "GlobalIntegral.h" +#include "GlobalNodes.h" #include "LocalIntegral.h" #include "IntegrandBase.h" #include "CoordinateMapping.h" @@ -363,10 +364,12 @@ void ASMu2D::extendRefinementDomain (IntSet& refineIndices, IntVec bndry0; std::vector bndry1(4); - for (int i = 0; i < 4; i++) + for (int i = 1; i <= 4; i++) { - bndry0.push_back(this->getCorner((i%2)*2-1, (i/2)*2-1, 1)); - this->getBoundaryNodes(1+i, bndry1[i], 1, 1, 0, true); + bndry0.push_back(GlobalNodes::getBoundaryNodes(*this->getBasis(ASM::REFINEMENT_BASIS), + 0, i, 0).front()); + bndry1[i-1] = GlobalNodes::getBoundaryNodes(*this->getBasis(ASM::REFINEMENT_BASIS), + 1, i, 0); } // Add refinement from neighbors @@ -377,7 +380,7 @@ void ASMu2D::extendRefinementDomain (IntSet& refineIndices, // Check if node is a corner node, // compute large extended domain (all directions) for (int edgeNode : bndry0) - if (edgeNode-1 == j) + if (edgeNode == j) { IntVec secondary = this->getOverlappingNodes(j); refineIndices.insert(secondary.begin(),secondary.end()); @@ -389,7 +392,7 @@ void ASMu2D::extendRefinementDomain (IntSet& refineIndices, // compute small extended domain (one direction) for (int edge = 0; edge < 4 && !done_with_this_node; edge++) for (int edgeNode : bndry1[edge]) - if (edgeNode-1 == j) + if (edgeNode == j) { IntVec secondary = this->getOverlappingNodes(j, edge/2+1); refineIndices.insert(secondary.begin(),secondary.end()); diff --git a/src/ASM/LR/ASMu3D.C b/src/ASM/LR/ASMu3D.C index 25b8e2aec..f815cbed9 100644 --- a/src/ASM/LR/ASMu3D.C +++ b/src/ASM/LR/ASMu3D.C @@ -21,6 +21,7 @@ #include "TimeDomain.h" #include "FiniteElement.h" #include "GlobalIntegral.h" +#include "GlobalNodes.h" #include "LocalIntegral.h" #include "IntegrandBase.h" #include "CoordinateMapping.h" @@ -2253,18 +2254,18 @@ void ASMu3D::extendRefinementDomain (IntSet& refineIndices, const int nface = 6; IntVec bndry0; - for (int K = -1; K < 2; K += 2) - for (int J = -1; J < 2; J += 2) - for (int I = -1; I < 2; I += 2) - bndry0.push_back(this->getCorner(I,J,K,1)); + const LR::LRSplineVolume* ref = this->getBasis(ASM::REFINEMENT_BASIS); + for (size_t c = 1; c <= 8; ++c) + bndry0.push_back(GlobalNodes::getBoundaryNodes(*ref, + 0, c, 0).front()); std::vector bndry1(nedge); for (int j = 1; j <= nedge; j++) - this->getBoundary1Nodes(j , bndry1[j-1], 1, 0, true); + bndry1[j-1] = GlobalNodes::getBoundaryNodes(*ref, 1, j, 0); std::vector bndry2(nface); for (int j = 1; j <= nface; j++) - this->getBoundaryNodes(j, bndry2[j-1], 1, 1, 0, true); + bndry2[j-1] = GlobalNodes::getBoundaryNodes(*ref, 2, j, 0); // Add refinement from neighbors for (int j : neighborIndices) @@ -2274,7 +2275,7 @@ void ASMu3D::extendRefinementDomain (IntSet& refineIndices, // Check if node is a corner node, // compute large extended domain (all directions) for (int edgeNode : bndry0) - if (edgeNode-1 == j) + if (edgeNode == j) { IntVec secondary = this->getOverlappingNodes(j); refineIndices.insert(secondary.begin(),secondary.end()); @@ -2287,7 +2288,7 @@ void ASMu3D::extendRefinementDomain (IntSet& refineIndices, int allowedDir; for (int edge = 0; edge < nedge && !done_with_this_node; edge++) for (int edgeNode : bndry1[edge]) - if (edgeNode-1 == j) + if (edgeNode == j) { if (edge < 4) allowedDir = 6; // bin(110), allowed to grow in v- and w-direction @@ -2305,7 +2306,7 @@ void ASMu3D::extendRefinementDomain (IntSet& refineIndices, // compute small extended domain (1 direction) for (int face = 0; face < nface && !done_with_this_node; face++) for (int edgeNode : bndry2[face]) - if (edgeNode-1 == j) + if (edgeNode == j) { allowedDir = 1 << face/2; IntVec secondary = this->getOverlappingNodes(j,allowedDir); diff --git a/src/SIM/AdaptiveSetup.C b/src/SIM/AdaptiveSetup.C index 4b4c1cdf3..d0e912fb9 100644 --- a/src/SIM/AdaptiveSetup.C +++ b/src/SIM/AdaptiveSetup.C @@ -21,7 +21,10 @@ #include "IFEM.h" #include "tinyxml2.h" #ifdef HAS_LRSPLINE +#include "ASMLRSpline.h" +#include "GlobalNodes.h" #include +#include #endif #include #include @@ -331,16 +334,21 @@ int AdaptiveSetup::calcRefinement (LR::RefineData& prm, int iStep, if (scheme == ISOTROPIC_FUNCTION) // use errors per function { - if (thePatch->getNoRefineNodes() == thePatch->getNoNodes(1)) - error.resize(model.getNoNodes(1),DblIdx(0.0,0)); - else if (model.getNoPatches() == 1) - error.resize(thePatch->getNoRefineNodes(),DblIdx(0.0,0)); - else - { - std::cerr <<" *** AdaptiveSetup::calcRefinement: Multi-patch refinement" - <<" is not available for mixed models."<< std::endl; - return -3; - } + size_t nNodes = 0; + if (model.getNoPatches() > 1) { +#ifdef HAS_LRSPLINE + std::vector refBasis; + for (const ASMbase* pch : model.getFEModel()) + refBasis.push_back(dynamic_cast(pch)->getRefinementBasis()); + + prm.MLGN = GlobalNodes::calcGlobalNodes(refBasis, model.getInterfaces()); + nNodes = *std::max_element(prm.MLGN.back().begin(), prm.MLGN.back().end()) + 1; +#endif + } else + nNodes = model.getPatch(1)->getNoRefineNodes(); + error.reserve(nNodes); + for (i = 0; i < nNodes; i++) + error.push_back(DblIdx(0.0,i)); for (i = 0; i < error.size(); i++) error[i].second = i; @@ -359,7 +367,7 @@ int AdaptiveSetup::calcRefinement (LR::RefineData& prm, int iStep, // Insert into global error array for (i = 0; i < locErr.size(); i++) if (model.getNoPatches() > 1) - error[patch->getNodeID(i+1)-1].first += locErr[i]; + error[prm.MLGN[patch->idx][i]].first += locErr[i]; else error[i].first += locErr[i]; } diff --git a/src/SIM/SIMinput.C b/src/SIM/SIMinput.C index bcf8f3eb9..2dd052e63 100644 --- a/src/SIM/SIMinput.C +++ b/src/SIM/SIMinput.C @@ -1604,9 +1604,16 @@ bool SIMinput::refine (const LR::RefineData& prm, Vectors& sol) // Extract local indices from the vector of global indices int locId; for (int k : prm.elements) - if ((locId = myModel[i]->getNodeIndex(k+1)) > 0) - if (refineIndices[i].insert(locId-1).second) + if (!prm.MLGN.empty()) { + const auto it = std::find(prm.MLGN[i].begin(), prm.MLGN[i].end(), k); + if (it != prm.MLGN[i].end()) + if (refineIndices[i].insert(std::distance(prm.MLGN[i].begin(), it)).second) + changed = true; + } else { + if ((locId = myModel[i]->getNodeIndex(k+1)) > 0) + if (refineIndices[i].insert(locId-1).second) changed = true; + } // Fetch boundary nodes covered (may need to pass this to other patches) IntVec bndry_nodes = pch[i]->getBoundaryCovered(refineIndices[i]); @@ -1627,13 +1634,27 @@ bool SIMinput::refine (const LR::RefineData& prm, Vectors& sol) for (int k : bndry_nodes) { // Check if this boundary node appears on other patches - int globId = myModel[i]->getNodeID(k+1); + int globId; + if (prm.MLGN.empty()) + globId = myModel[i]->getNodeID(k+1); + else + globId = prm.MLGN[i][k]; for (size_t j = 0; j < myModel.size(); j++) - if (j != i && (locId = myModel[j]->getNodeIndex(globId)) > 0) - if (conformingIndices[j].insert(locId-1).second) { - changed = true; - conformingIndices[i].insert(k); + if (j != i) { + if (prm.MLGN.empty()) { + if ((locId = myModel[j]->getNodeIndex(globId)) > 0) { + conformingIndices[j].insert(locId-1); + conformingIndices[i].insert(k); + } + } else { + const auto it = std::find(prm.MLGN[j].begin(), + prm.MLGN[j].end(), globId); + if (it != prm.MLGN[j].end()) { + conformingIndices[j].insert(std::distance(prm.MLGN[j].begin(), it)); + conformingIndices[i].insert(k); + } } + } } } diff --git a/src/SIM/SIMinput.h b/src/SIM/SIMinput.h index 392b18a66..862b5cea4 100644 --- a/src/SIM/SIMinput.h +++ b/src/SIM/SIMinput.h @@ -208,6 +208,10 @@ class SIMinput : public SIMbase coordCheck); } + //! \brief Obtain a const reference to model topology. + const std::vector& getInterfaces() const + { return myInterfaces; } + protected: //! \brief Helper method returning a stream for patch geometry input. //! \param[in] tag The XML-tag containing the patch geometry definition