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 b966e25a8..9166cb291 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/ASMLRSpline.h b/src/ASM/LR/ASMLRSpline.h index e6c629a63..62124477c 100644 --- a/src/ASM/LR/ASMLRSpline.h +++ b/src/ASM/LR/ASMLRSpline.h @@ -108,6 +108,9 @@ class ASMLRSpline : public ASMbase, public ASMunstruct //! \brief Returns parameter values and node numbers of the domain corners. virtual bool getParameterDomain(Real2DMat&, IntVec*) const; + //! \brief Returns a const pointer to refinement basis. + const LR::LRSpline* getRefinementBasis() const { return refB.get(); } + //! \brief Refines the mesh adaptively. //! \param[in] prm Input data used to control the mesh refinement //! \param sol Control point results values that are transferred to new mesh diff --git a/src/ASM/LR/ASMu2D.C b/src/ASM/LR/ASMu2D.C index 9291ab336..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()); @@ -2464,9 +2467,6 @@ Field* ASMu2D::getProjectedField (const Vector& coefs) const Fields* ASMu2D::getProjectedFields (const Vector& coefs, size_t) const { - if (!this->separateProjectionBasis()) - return nullptr; - size_t ncmp = coefs.size() / this->getNoProjectionNodes(); if (ncmp*this->getNoProjectionNodes() == coefs.size()) return new LRSplineFields2D(this->getBasis(ASM::PROJECTION_BASIS), diff --git a/src/ASM/LR/ASMu3D.C b/src/ASM/LR/ASMu3D.C index 90e0d8f35..dd1d2c699 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" @@ -2004,9 +2005,6 @@ Field* ASMu3D::getProjectedField (const Vector& coefs) const Fields* ASMu3D::getProjectedFields (const Vector& coefs, size_t) const { - if (!this->separateProjectionBasis()) - return nullptr; - size_t ncmp = coefs.size() / this->getNoProjectionNodes(); if (ncmp*this->getNoProjectionNodes() == coefs.size()) return new LRSplineFields3D(this->getBasis(ASM::PROJECTION_BASIS),coefs,ncmp); @@ -2254,18 +2252,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) @@ -2275,7 +2273,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()); @@ -2288,7 +2286,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 @@ -2306,7 +2304,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/ASM/LR/GlobalNodes.C b/src/ASM/LR/GlobalNodes.C new file mode 100644 index 000000000..349e5705e --- /dev/null +++ b/src/ASM/LR/GlobalNodes.C @@ -0,0 +1,186 @@ +// $Id$ +//============================================================================== +//! +//! \file GlobalNodes.C +//! +//! \date Mar 13 2018 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Simple global node establishment for unstructured FE models. +//! +//============================================================================== + +#include "GlobalNodes.h" + +#include "ASMLRSpline.h" +#include "Utilities.h" + +#include + +namespace { + +/*! + \brief Class for ordering interfaces in processing order. +*/ + +class InterfaceOrder { +public: + //! \brief Comparison operator for interfaces + //! \param A First interface + //! \param B Second interface + bool operator()(const ASM::Interface& A, const ASM::Interface& B) const + { + if (A.master != B.master) + return A.master < B.master; + + if (A.slave != B.slave) + return A.slave < B.slave; + + if (A.dim != B.dim) + return A.dim < B.dim; + + return A.midx < B.midx; + } +}; + +} + + +GlobalNodes::IntVec +GlobalNodes::getBoundaryNodes (const LR::LRSpline& lr, + int dim, int lidx, int orient) +{ + LR::parameterEdge edge{}; + if (dim == 0) { + if (lr.nVariate() == 2) { + switch (lidx) { + case 1: edge = LR::WEST | LR::SOUTH; break; + case 2: edge = LR::EAST | LR::SOUTH; break; + case 3: edge = LR::WEST | LR::NORTH; break; + case 4: edge = LR::EAST | LR::NORTH; break; + } + } else { + switch (lidx) { + case 1: edge = LR::WEST | LR::SOUTH | LR::BOTTOM; break; + case 2: edge = LR::EAST | LR::SOUTH | LR::BOTTOM; break; + case 3: edge = LR::WEST | LR::NORTH | LR::BOTTOM; break; + case 4: edge = LR::EAST | LR::NORTH | LR::BOTTOM; break; + case 5: edge = LR::WEST | LR::SOUTH | LR::TOP; break; + case 6: edge = LR::EAST | LR::SOUTH | LR::TOP; break; + case 7: edge = LR::WEST | LR::NORTH | LR::TOP; break; + case 8: edge = LR::EAST | LR::NORTH | LR::TOP; break; + } + } + } else if (dim == 1) { + if (lr.nVariate() == 2) { + switch (lidx) { + case 1: edge = LR::WEST; break; + case 2: edge = LR::EAST; break; + case 3: edge = LR::SOUTH; break; + case 4: edge = LR::NORTH; break; + default: break; + } + } else { + switch (lidx) { + case 1: edge = LR::BOTTOM | LR::SOUTH; break; + case 2: edge = LR::BOTTOM | LR::NORTH; break; + case 3: edge = LR::TOP | LR::SOUTH; break; + case 4: edge = LR::TOP | LR::NORTH; break; + case 5: edge = LR::BOTTOM | LR::WEST; break; + case 6: edge = LR::BOTTOM | LR::EAST; break; + case 7: edge = LR::TOP | LR::WEST; break; + case 8: edge = LR::TOP | LR::WEST; break; + case 9: edge = LR::SOUTH | LR::WEST; break; + case 10: edge = LR::SOUTH | LR::EAST; break; + case 11: edge = LR::NORTH | LR::WEST; break; + case 12: edge = LR::NORTH | LR::EAST; break; + } + } + } else if (dim == 2) { + switch (lidx) { + case 1: edge = LR::WEST; break; + case 2: edge = LR::EAST; break; + case 3: edge = LR::SOUTH; break; + case 4: edge = LR::NORTH; break; + case 5: edge = LR::BOTTOM; break; + case 6: edge = LR::TOP; break; + } + } + + std::vector edgeFunctions; + lr.getEdgeFunctions(edgeFunctions, edge); + + if (dim == 1) { + if (lr.nVariate() == 2) { + int v = (lidx == 1 || lidx == 2) ? 0 : 1; + int u = 1-v; + ASMLRSpline::Sort(u, v, orient, edgeFunctions); + } else { + int dir = (lidx-1)/4; + int u = dir == 0; + int v = 1 + (dir != 2); + ASMLRSpline::Sort(u, v, orient, edgeFunctions); + } + } else if (dim == 2) { + int dir = (lidx-1)/2; + int u = dir == 0; + int v = 1 + (dir != 2); + ASMLRSpline::Sort(u, v, orient, edgeFunctions); + } + + GlobalNodes::IntVec lNodes; + lNodes.reserve(edgeFunctions.size()); + for (const LR::Basisfunction* func : edgeFunctions) + lNodes.push_back(func->getId()); + + return lNodes; +} + + +std::vector +GlobalNodes::calcGlobalNodes (const GlobalNodes::LRSplineVec& pchs, + const GlobalNodes::InterfaceVec& interfaces) +{ + // count total number of nodes + size_t nNodes = 0; + std::vector result(pchs.size()); + auto it = result.begin(); + for (const auto& pch : pchs) { + it->resize(pch->nBasisFunctions()); + std::iota(it->begin(), it->end(), nNodes); + nNodes += pch->nBasisFunctions(); + ++it; + } + + // remap common nodes + std::set ifset; + for (const ASM::Interface& it : interfaces) + ifset.insert(it); + for (size_t i = 0; i < pchs.size(); ++i) { + std::map old2new; + for (const ASM::Interface& it : ifset) { + if (it.master != (int)i+1) + continue; + + IntVec mNodes = getBoundaryNodes(*pchs[i], it.dim, it.midx, 0); + IntVec sNodes = getBoundaryNodes(*pchs[it.slave-1], it.dim, it.sidx, it.orient); + for (size_t n = 0; n < mNodes.size(); ++n) + old2new[result[it.slave-1][sNodes[n]]] = result[i][mNodes[n]]; + } + + // renumber + for (size_t j = i; j < pchs.size(); ++j) + for (int& it : result[j]) + utl::renumber(it, old2new, false); + + // compress + int maxNode = *std::max_element(result[i].begin(), result[i].end()); + for (size_t j = i+1; j < pchs.size(); ++j) + for (int& n : result[j]) + if (n > maxNode) + n = ++maxNode; + } + + return result; +} diff --git a/src/ASM/LR/GlobalNodes.h b/src/ASM/LR/GlobalNodes.h new file mode 100644 index 000000000..3e7a38461 --- /dev/null +++ b/src/ASM/LR/GlobalNodes.h @@ -0,0 +1,50 @@ +// $Id$ +//============================================================================== +//! +//! \file GlobalNodes.h +//! +//! \date Mar 13 2018 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Simple global node establishment for unstructured FE models. +//! +//============================================================================== + +#ifndef _GLOBAL_NODES_H_ +#define _GLOBAL_NODES_H_ + +#include "Interface.h" +#include + +#include + + +/*! + \brief Class establishing global node numbers for unstructed FE models. +*/ + +class GlobalNodes +{ +public: + using IntVec = std::vector; //!< Convenience typedef + using LRSplineVec = std::vector; //!< Convenience typedef + using InterfaceVec = std::vector; //!< Convenience typedef + + //! \brief Extract local boundary nodes for a LR spline. + //! \param lr The LR spline to extract boundary nodes for + //! \param dim The dimension of the boundary to extract + //! \param lidx The local index of the boundary to extract + //! \param orient Orientation of nodes on boundary + static IntVec getBoundaryNodes(const LR::LRSpline& lr, + int dim, int lidx, int orient); + + + //! \brief Calculate global node numbers for a FE model. + //! \param pchs The spline patches in the model + //! \param interfaces The topological connections for the spline patches + static std::vector calcGlobalNodes(const LRSplineVec& pchs, + const InterfaceVec& interfaces); +}; + +#endif diff --git a/src/ASM/LR/Test/ASMuCube.h b/src/ASM/LR/Test/ASMuCube.h index 4091abed7..cfe4d32b2 100644 --- a/src/ASM/LR/Test/ASMuCube.h +++ b/src/ASM/LR/Test/ASMuCube.h @@ -43,14 +43,30 @@ class ASMuCube : public ASMu3D "# Elements:\n" "0 [3] : (0, 0, 0) x (1, 1, 1) {0, 1, 2, 3, 4, 5, 6, 7}\n"; - explicit ASMuCube(unsigned char n_f = 3) : ASMu3D(n_f) + explicit ASMuCube(unsigned char n_f = 3, + double xshift = 0.0, + double yshift = 0.0, + double zshift = 0.0) + : ASMu3D(n_f) { + auto&& addLine = [](double x, double y, double z) + { + return std::to_string(x) + " " + + std::to_string(y) + " " + + std::to_string(z) + '\n'; + }; std::stringstream geo("700 1 0 0\n3 0\n" "2 2\n0 0 1 1\n" "2 2\n0 0 1 1\n" - "2 2\n0 0 1 1\n" - "0 0 0\n1 0 0\n0 1 0\n1 1 0\n" - "0 0 1\n1 0 1\n0 1 1\n1 1 1\n"); + "2 2\n0 0 1 1\n" + + addLine(xshift, yshift, zshift) + + addLine(xshift + 1.0, yshift, zshift) + + addLine(xshift, yshift + 1.0, zshift) + + addLine(xshift + 1.0, yshift + 1.0, zshift) + + addLine(xshift, yshift, zshift + 1.0) + + addLine(xshift + 1.0, yshift, zshift + 1.0) + + addLine(xshift, yshift + 1.0, zshift + 1.0) + + addLine(xshift + 1.0, yshift + 1.0, zshift + 1.0)); this->read(geo); } virtual ~ASMuCube() {} diff --git a/src/ASM/LR/Test/ASMuSquare.h b/src/ASM/LR/Test/ASMuSquare.h index e6a3bd014..f03c9a4d7 100644 --- a/src/ASM/LR/Test/ASMuSquare.h +++ b/src/ASM/LR/Test/ASMuSquare.h @@ -34,9 +34,21 @@ class ASMuSquare : public ASMu2D "[0, 1] x 1 (2)\n" "# Elements:\n" "0 [2] : (0, 0) x (1, 1) {0, 1, 2, 3}\n"; - explicit ASMuSquare(unsigned char n_f = 2) : ASMu2D(2,n_f) + + explicit ASMuSquare(unsigned char n_f = 2, + double xshift = 0.0, + double yshift = 0.0) + : ASMu2D(2,n_f) { - std::stringstream geo("200 1 0 0\n2 0\n2 2\n0 0 1 1\n2 2\n0 0 1 1\n0 0\n1 0\n0 1\n1 1\n"); + auto&& addLine = [](double x, double y) + { + return std::to_string(x) + " " + std::to_string(y) + '\n'; + }; + std::stringstream geo("200 1 0 0\n2 0\n2 2\n0 0 1 1\n2 2\n0 0 1 1\n" + + addLine(xshift, yshift) + + addLine(xshift+1.0, yshift) + + addLine(xshift, yshift+1.0) + + addLine(xshift+1.0, yshift+1.0)); this->read(geo); } virtual ~ASMuSquare() {} diff --git a/src/ASM/LR/Test/TestGlobalNodes.C b/src/ASM/LR/Test/TestGlobalNodes.C new file mode 100644 index 000000000..e705fd3a3 --- /dev/null +++ b/src/ASM/LR/Test/TestGlobalNodes.C @@ -0,0 +1,74 @@ +//============================================================================== +//! +//! \file TestGlobalNodes.C +//! +//! \date Jun 21 2024 +//! +//! \author Arne Morten Kvarving / SINTEF +//! +//! \brief Tests for simple global node establishment for unstructured FE models. +//! +//============================================================================== + +#include "GlobalNodes.h" +#include "ASMuCube.h" +#include "ASMuSquare.h" + +#include "LRSpline/LRSplineVolume.h" + +#include "gtest/gtest.h" + + +TEST(TestGlobalNodes, 2D) +{ + ASMuSquare pch1; + pch1.generateFEMTopology(); + ASMuSquare pch2(2, 1.0, 0.0); + pch2.generateFEMTopology(); + ASMuSquare pch3(2, 1.0, 1.0); + pch3.generateFEMTopology(); + + std::vector splines{pch1.getBasis(), + pch2.getBasis(), + pch3.getBasis()}; + std::vector ifs; + ifs.push_back(ASM::Interface{1, 2, 2, 1, 0, 1, 1, 0}); + ifs.push_back(ASM::Interface{2, 3, 4, 3, 0, 1, 1, 0}); + + auto nodes = GlobalNodes::calcGlobalNodes(splines, ifs); + const auto ref = std::vector{ + GlobalNodes::IntVec{0, 1, 2, 3}, + GlobalNodes::IntVec{1, 4, 3, 5}, + GlobalNodes::IntVec{3, 5, 6, 7}, + }; + + EXPECT_EQ(nodes, ref); +} + + +TEST(TestGlobalNodes, 3D) +{ + ASMuCube pch1; + pch1.generateFEMTopology(); + ASMuCube pch2(2, 1.0, 0.0, 0.0); + pch2.generateFEMTopology(); + ASMuCube pch3(2, 1.0, 1.0, 1.0); + pch3.generateFEMTopology(); + + std::vector splines{pch1.getBasis(), + pch2.getBasis(), + pch3.getBasis()}; + std::vector ifs; + ifs.push_back(ASM::Interface{1, 2, 2, 1, 0, 2, 1, 0}); + ifs.push_back(ASM::Interface{2, 3, 6, 5, 0, 2, 1, 0}); + + auto nodes = GlobalNodes::calcGlobalNodes(splines, ifs); + + const auto ref = std::vector{ + GlobalNodes::IntVec{0, 1, 2, 3, 4, 5, 6, 7}, + GlobalNodes::IntVec{1, 8, 3, 9, 5, 10, 7, 11}, + GlobalNodes::IntVec{5, 10, 7, 11, 12, 13, 14, 15}, + }; + + EXPECT_EQ(nodes, ref); +} diff --git a/src/SIM/AdaptiveSetup.C b/src/SIM/AdaptiveSetup.C index 912862d46..3627662af 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/SIMbase.C b/src/SIM/SIMbase.C index 56a42d2e5..a3ef6ef64 100644 --- a/src/SIM/SIMbase.C +++ b/src/SIM/SIMbase.C @@ -2119,7 +2119,7 @@ bool SIMbase::project (Matrix& ssol, const Vector& psol, ssol.resize(myProblem->getNoFields(2),ngNodes); ssol.fillBlock(values, 1, ofs+1); - ofs += myModel[i]->getNoProjectionNodes()*myProblem->getNoFields(2); + ofs += myModel[i]->getNoProjectionNodes(); } else { size_t nComps = values.rows(); size_t nNodes = values.cols(); @@ -2369,6 +2369,11 @@ bool SIMbase::evalSecondarySolution (Matrix& field, int pindx) const bool SIMbase::fieldProjections () const { + // Forced for mixed multi-patch models as the nodal averaging + // in SIMbase::project() fails. + if (myModel.size() > 1 && myModel[0]->getNoBasis() > 1) + return true; + for (const ASMbase* pch : myModel) if (pch->separateProjectionBasis()) return true; diff --git a/src/SIM/SIMinput.C b/src/SIM/SIMinput.C index 9c50d8ad9..017fcae40 100644 --- a/src/SIM/SIMinput.C +++ b/src/SIM/SIMinput.C @@ -1606,9 +1606,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) - changed = true; + 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]); @@ -1629,13 +1636,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