From b7f97c498ce4be5b20eaa784e1ed229ecd98922f Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 14 Mar 2018 15:59:04 +0100 Subject: [PATCH] added: class for simple global node number establishment for LR splines use this to enable adaptive refinement for mixed multi-patch models. in particular, the refinements basis differ from the finite element bases for some formulations, so a simple node numbering scheme for the refinement basis is required to perform the multi-patch refinements. --- CMakeLists.txt | 3 + src/ASM/ASMunstruct.h | 4 + src/ASM/GlobalNodes.C | 173 +++++++++++++++++++++++++++++++++++++++ src/ASM/GlobalNodes.h | 50 +++++++++++ src/ASM/LR/ASMLRSpline.C | 15 ++-- src/ASM/LR/ASMu2D.C | 13 +-- src/ASM/LR/ASMu2Dmx.C | 6 ++ src/ASM/LR/ASMu2Dmx.h | 3 + src/ASM/LR/ASMu3D.C | 59 +++++++------ src/ASM/LR/ASMu3D.h | 6 +- src/ASM/LR/ASMu3Dmx.C | 29 ++++++- src/ASM/LR/ASMu3Dmx.h | 18 ++++ src/SIM/AdaptiveSIM.C | 33 +++++--- src/SIM/SIMbase.C | 2 +- src/SIM/SIMinput.C | 38 +++++++-- src/SIM/SIMinput.h | 3 + 16 files changed, 387 insertions(+), 68 deletions(-) create mode 100644 src/ASM/GlobalNodes.C create mode 100644 src/ASM/GlobalNodes.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 4fbc4cc28..33d7c1b9b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -104,6 +104,9 @@ endif() if(LRSPLINE_FOUND OR LRSpline_FOUND) file(GLOB LR_SRCS ${PROJECT_SOURCE_DIR}/src/ASM/LR/*.C) list(APPEND IFEM_SRCS ${LR_SRCS}) +else() + list(REMOVE_ITEM IFEM_SRCS + ${PROJECT_SOURCE_DIR}/src/ASM/GlobalNodes.C) endif() list(APPEND CHECK_SOURCES ${IFEM_SRCS}) add_library(IFEM ${IFEM_SRCS} ${TINYXML_SRCS}) diff --git a/src/ASM/ASMunstruct.h b/src/ASM/ASMunstruct.h index 5ef29e867..474015881 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) {} @@ -214,6 +215,9 @@ class ASMunstruct : public ASMbase //! \param[in] refTol Mesh refinement threshold virtual bool refine(const RealFunc& refC, double refTol) = 0; + //! \brief Obtain the refinement basis. + virtual const LR::LRSpline* getRefinementBasis() const { return geo; } + protected: //! \brief Refines the mesh adaptively. //! \param[in] prm Input data used to control the mesh refinement diff --git a/src/ASM/GlobalNodes.C b/src/ASM/GlobalNodes.C new file mode 100644 index 000000000..bdd48d5ca --- /dev/null +++ b/src/ASM/GlobalNodes.C @@ -0,0 +1,173 @@ +// $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 "ASMunstruct.h" +#include "Utilities.h" + + +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; + ASMunstruct::Sort(u, v, orient, edgeFunctions); + } else { + int dir = (lidx-1)/4; + int u = dir == 0; + int v = 1 + (dir != 2); + ASMunstruct::Sort(u, v, orient, edgeFunctions); + } + } else if (dim == 2) { + int dir = (lidx-1)/2; + int u = dir == 0; + int v = 1 + (dir != 2); + ASMunstruct::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; +} + + +class InterfaceOrder { +public: + 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; + } +}; + + +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 LR::LRSpline* pch : pchs) { + it->resize(pch->nBasisFunctions()); + std::iota(it->begin(), it->end(), nNodes); + nNodes += pch->nBasisFunctions(); + ++it; + } + + // remap common nodes + InterfaceOrder ifOrder; + std::set ifset(ifOrder); + 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/GlobalNodes.h b/src/ASM/GlobalNodes.h new file mode 100644 index 000000000..27a1b3d21 --- /dev/null +++ b/src/ASM/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: + typedef std::vector IntVec; //!< Convenience typedef + typedef std::vector LRSplineVec; //!< Convenience typedef + typedef std::vector InterfaceVec; //!< 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/ASMLRSpline.C b/src/ASM/LR/ASMLRSpline.C index f9f3676cb..7c87e827a 100644 --- a/src/ASM/LR/ASMLRSpline.C +++ b/src/ASM/LR/ASMLRSpline.C @@ -15,6 +15,7 @@ #include "IFEM.h" #include "LRSpline/LRSplineSurface.h" #include "LRSpline/Basisfunction.h" +#include "GlobalNodes.h" #include "Profiler.h" #include "ThreadGroups.h" #include @@ -334,15 +335,16 @@ void ASMunstruct::getFunctionsForElements (IntSet& functions, IntVec ASMunstruct::getBoundaryNodesCovered (const IntSet& nodes) const { IntSet result; + const LR::LRSpline* refBasis = this->getRefinementBasis(); int numbEdges = (this->getNoParamDim() == 2) ? 4 : 6; for (int edge = 1; edge <= numbEdges; edge++) { - IntVec oneBoundary; - this->getBoundaryNodes(edge, oneBoundary, 1, 1, 0, true); // this returns a 1-indexed list + IntVec bnd = GlobalNodes::getBoundaryNodes(*refBasis, + this->getNoParamDim()-1, edge, 0); for (const int i : nodes) - for (const int j : oneBoundary) - if (geo->getBasisfunction(i)->contains(*geo->getBasisfunction(j-1))) - result.insert(j-1); + for (const int j : bnd) + if (refBasis->getBasisfunction(i)->contains(*refBasis->getBasisfunction(j))) + result.insert(j); } return IntVec(result.begin(), result.end()); @@ -352,9 +354,10 @@ IntVec ASMunstruct::getBoundaryNodesCovered (const IntSet& nodes) const IntVec ASMunstruct::getOverlappingNodes (const IntSet& nodes, int dir) const { IntSet result; + const LR::LRSpline* refBasis = this->getRefinementBasis(); for (const int i : nodes) { - LR::Basisfunction *b = geo->getBasisfunction(i); + const LR::Basisfunction *b = refBasis->getBasisfunction(i); for (auto el : b->support()) // for all elements where *b has support for (auto basis : el->support()) // for all functions on this element { diff --git a/src/ASM/LR/ASMu2D.C b/src/ASM/LR/ASMu2D.C index 7115323ed..1202a000e 100644 --- a/src/ASM/LR/ASMu2D.C +++ b/src/ASM/LR/ASMu2D.C @@ -22,6 +22,7 @@ #include "TimeDomain.h" #include "FiniteElement.h" #include "GlobalIntegral.h" +#include "GlobalNodes.h" #include "LocalIntegral.h" #include "IntegrandBase.h" #include "CoordinateMapping.h" @@ -355,10 +356,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->getRefinementBasis(), + 0, i, 0).front()); + bndry1[i-1] = GlobalNodes::getBoundaryNodes(*this->getRefinementBasis(), + 1, i, 0); } // Add refinement from neighbors @@ -369,7 +372,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()); @@ -381,7 +384,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); refineIndices.insert(secondary.begin(),secondary.end()); diff --git a/src/ASM/LR/ASMu2Dmx.C b/src/ASM/LR/ASMu2Dmx.C index d7e86276f..80970e031 100644 --- a/src/ASM/LR/ASMu2Dmx.C +++ b/src/ASM/LR/ASMu2Dmx.C @@ -1200,3 +1200,9 @@ bool ASMu2Dmx::evalProjSolution (Matrix& sField, const Vector& locSol, return true; } + + +const LR::LRSpline* ASMu2Dmx::getRefinementBasis() const +{ + return refBasis.get(); +} diff --git a/src/ASM/LR/ASMu2Dmx.h b/src/ASM/LR/ASMu2Dmx.h index a3073d91d..3c2f2d566 100644 --- a/src/ASM/LR/ASMu2Dmx.h +++ b/src/ASM/LR/ASMu2Dmx.h @@ -222,6 +222,9 @@ class ASMu2Dmx : public ASMu2D, private ASMmxBase virtual void remapErrors(RealArray& errors, const RealArray& origErr, bool elemErrors) const; + //! \brief Obtain the refinement basis. + virtual const LR::LRSpline* getRefinementBasis() const; + protected: //! \brief Assembles L2-projection matrices for the secondary solution. //! \param[out] A Left-hand-side matrix diff --git a/src/ASM/LR/ASMu3D.C b/src/ASM/LR/ASMu3D.C index a608f2d8c..3ce855b55 100644 --- a/src/ASM/LR/ASMu3D.C +++ b/src/ASM/LR/ASMu3D.C @@ -22,6 +22,7 @@ #include "TimeDomain.h" #include "FiniteElement.h" #include "GlobalIntegral.h" +#include "GlobalNodes.h" #include "LagrangeInterpolator.h" #include "LocalIntegral.h" #include "IntegrandBase.h" @@ -276,7 +277,7 @@ bool ASMu3D::connectPatch (int face, ASM3D& neighbor, int nface, bool ASMu3D::connectBasis (int face, ASMu3D& neighbor, int nface, int norient, int basis, int slave, int master, - bool coordCheck, int thick) + bool coordCheck, int) { if (this->isShared() && neighbor.isShared()) return true; @@ -305,28 +306,23 @@ bool ASMu3D::connectBasis (int face, ASMu3D& neighbor, int nface, int norient, } const double xtol = 1.0e-4; - for (size_t node = 0; node < masterNodes.size(); ++node) - for (int t = 0; t < thick; ++t) + for (size_t node = 0; node < masterNodes.size(); ++node) { + int node2 = masterNodes[node]; + int slave = slaveNodes[node]; + + if (!coordCheck) + ASMbase::collapseNodes(neighbor,node2,*this,slave); + else if (neighbor.getCoord(node2).equal(this->getCoord(slave),xtol)) + ASMbase::collapseNodes(neighbor,node2,*this,slave); + else { - // TODO: Hey, this looks suspicious. If thick > 1, this will clearly read - // outside the master/slaveNodes arrays. If thick > 1 is not supported - // for LR, please remove it alltogether in this class to avoid confusion. - int node2 = masterNodes[node*thick+t]; - int slave = slaveNodes[node*thick+t]; - - if (!coordCheck) - ASMbase::collapseNodes(neighbor,node2,*this,slave); - else if (neighbor.getCoord(node2).equal(this->getCoord(slave),xtol)) - ASMbase::collapseNodes(neighbor,node2,*this,slave); - else - { - std::cerr <<" *** ASMu3D::connectPatch: Non-matching nodes " - << node2 <<": "<< neighbor.getCoord(node2) - <<"\n and " - << slave <<": "<< this->getCoord(slave) << std::endl; - return false; - } + std::cerr <<" *** ASMu3D::connectPatch: Non-matching nodes " + << node2 <<": "<< neighbor.getCoord(node2) + <<"\n and " + << slave <<": "<< this->getCoord(slave) << std::endl; + return false; } + } return true; } @@ -2478,18 +2474,19 @@ 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)); + for (size_t c = 1; c <= 8; ++c) + bndry0.push_back(GlobalNodes::getBoundaryNodes(*this->getRefinementBasis(), + 0, c, 0).front()); std::vector bndry1; for (int j = 1; j <= nedge; j++) - bndry1.push_back(this->getEdge(j, true, 1, 0)); + bndry1.push_back(GlobalNodes::getBoundaryNodes(*this->getRefinementBasis(), + 1, j, 0)); - std::vector bndry2(nface); + std::vector bndry2; for (int j = 1; j <= nface; j++) - this->getBoundaryNodes(j, bndry2[j-1], 1, 1, 0, true); + bndry2.push_back(GlobalNodes::getBoundaryNodes(*this->getRefinementBasis(), + 2, j, 0)); // Add refinement from neighbors for (int j : neighborIndices) @@ -2499,7 +2496,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()); @@ -2512,7 +2509,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 @@ -2530,7 +2527,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/ASMu3D.h b/src/ASM/LR/ASMu3D.h index 2f7935956..1b461640e 100644 --- a/src/ASM/LR/ASMu3D.h +++ b/src/ASM/LR/ASMu3D.h @@ -228,7 +228,6 @@ class ASMu3D : public ASMunstruct, public ASM3D //! \param[in] nface Local face index of neighbor patch, in range [1,6] //! \param[in] norient Relative face orientation flag (see below) //! \param[in] coordCheck False to disable coordinate checks (periodic connections) - //! \param[in] thick Thickness of connection //! //! \details The face orientation flag \a norient must be in range [0,7]. //! When interpreted as a binary number, its 3 digits are decoded as follows: @@ -236,7 +235,7 @@ class ASMu3D : public ASMunstruct, public ASM3D //! - middle digit = 1: Parameter \a u in neighbor patch face is reversed //! - right digit = 1: Parameter \a v in neighbor patch face is reversed virtual bool connectPatch(int face, ASM3D& neighbor, int nface, int norient, - int = 0, bool coordCheck = true, int thick = 1); + int = 0, bool coordCheck = true, int = 1); /* More multipatch stuff, maybe later... //! \brief Makes two opposite boundary faces periodic. @@ -463,10 +462,9 @@ class ASMu3D : public ASMunstruct, public ASM3D //! \param[in] slave 0-based index of the first slave node in this basis //! \param[in] master 0-based index of the first master node in this basis //! \param[in] coordCheck False to disable coordinate checks - //! \param[in] thick Thickness of connection bool connectBasis(int face, ASMu3D& neighbor, int nface, int norient, int basis = 1, int slave = 0, int master = 0, - bool coordCheck = true, int thick = 1); + bool coordCheck = true, int = 1); //! \brief Extracts parameter values of the Gauss points in one direction. //! \param[out] uGP Parameter values in given direction for all points diff --git a/src/ASM/LR/ASMu3Dmx.C b/src/ASM/LR/ASMu3Dmx.C index 15b975e2d..55ac45b8f 100644 --- a/src/ASM/LR/ASMu3Dmx.C +++ b/src/ASM/LR/ASMu3Dmx.C @@ -174,6 +174,9 @@ bool ASMu3Dmx::getSolution (Matrix& sField, const Vector& locSol, bool ASMu3Dmx::generateFEMTopology () { + if (!myMLGN.empty()) + return true; + if (m_basis.empty()) { auto vec = ASMmxBase::establishBases(tensorspline, ASMmxBase::Type); m_basis.resize(vec.size()); @@ -371,7 +374,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand, } int iEl = el->getId(); MxFiniteElement fe(elem_sizes); - fe.iel = iEl+1; + fe.iel = MLGE[iEl]; std::vector dNxdu(m_basis.size()); Matrix Xnod, Jac; @@ -508,7 +511,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand, else evaluateBasis(fe, dNxdu[b], bezierExtract[b][els[b]-1], B, b+1) ; } - fe.iel = iEl+1; + fe.iel = MLGE[iEl]; // Compute Jacobian inverse of coordinate mapping and derivatives fe.detJxW = utl::Jacobian(Jac,fe.grad(geoBasis),Xnod,dNxdu[geoBasis-1]); @@ -1086,3 +1089,25 @@ bool ASMu3Dmx::evalProjSolution (Matrix& sField, const Vector& locSol, return true; } + + +bool ASMu3Dmx::connectPatch (int face, ASM3D& neighbor, int nface, int norient, + int basis, bool coordCheck, int) +{ + ASMu3Dmx* neighMx = dynamic_cast(&neighbor); + if (!neighMx) return false; + + for (size_t i = 1; i <= m_basis.size(); ++i) + if (basis == 0 || i == (size_t)basis) + if (!this->connectBasis(face,*neighMx,nface,norient,i,0,0,coordCheck)) + return false; + + this->addNeighbor(neighMx); + return true; +} + + +const LR::LRSpline* ASMu3Dmx::getRefinementBasis() const +{ + return refBasis.get(); +} diff --git a/src/ASM/LR/ASMu3Dmx.h b/src/ASM/LR/ASMu3Dmx.h index 3c2345d99..58b272cb6 100644 --- a/src/ASM/LR/ASMu3Dmx.h +++ b/src/ASM/LR/ASMu3Dmx.h @@ -192,6 +192,24 @@ class ASMu3Dmx : public ASMu3D, private ASMmxBase virtual void remapErrors(RealArray& errors, const RealArray& origErr, bool elemErrors) const; + //! \brief Connects all matching nodes on two adjacent boundary faces. + //! \param[in] face Local face index of this patch, in range [1,6] + //! \param neighbor The neighbor patch + //! \param[in] nface Local face index of neighbor patch, in range [1,6] + //! \param[in] norient Relative face orientation flag (see below) + //! \param[in] coordCheck False to disable coordinate checks (periodic connections) + //! + //! \details The face orientation flag \a norient must be in range [0,7]. + //! When interpreted as a binary number, its 3 digits are decoded as follows: + //! - left digit = 1: The u and v parameters of the neighbor face are swapped + //! - middle digit = 1: Parameter \a u in neighbor patch face is reversed + //! - right digit = 1: Parameter \a v in neighbor patch face is reversed + virtual bool connectPatch(int face, ASM3D& neighbor, int nface, int norient, + int = 0, bool coordCheck = true, int = 1); + + //! \brief Obtain the refinement basis. + virtual const LR::LRSpline* getRefinementBasis() const; + protected: //! \brief Assembles L2-projection matrices for the secondary solution. //! \param[out] A Left-hand-side matrix diff --git a/src/SIM/AdaptiveSIM.C b/src/SIM/AdaptiveSIM.C index 04d669a67..e3ce6dbc0 100644 --- a/src/SIM/AdaptiveSIM.C +++ b/src/SIM/AdaptiveSIM.C @@ -16,6 +16,9 @@ #include "SIMenums.h" #include "ASMunstruct.h" #include "ASMmxBase.h" +#ifdef HAS_LRSPLINE +#include "GlobalNodes.h" +#endif #include "IntegrandBase.h" #include "Utilities.h" #include "IFEM.h" @@ -348,13 +351,19 @@ bool AdaptiveSIM::adaptMesh (int iStep) { if (model.getFEModel()[0]->getNoRefineNodes() != model.getFEModel()[0]->getNoNodes(1)) { + size_t nNodes = model.getFEModel()[0]->getNoRefineNodes(); if (model.getNoPatches() > 1) { - std::cerr <<" *** AdaptiveSIM::adaptMesh: Multi-patch refinement" - <<" is not available for mixed models."<< std::endl; - return false; + std::vector refBasis; + for (ASMbase* pch : model.getFEModel()) + refBasis.push_back(dynamic_cast(pch)->getRefinementBasis()); + +#ifdef HAS_LRSPLINE + prm.MLGN = GlobalNodes::calcGlobalNodes(refBasis, model.getInterfaces()); + nNodes = *std::max_element(prm.MLGN.back().begin(), prm.MLGN.back().end()) + 1; +#endif } - errors.reserve(model.getFEModel()[0]->getNoRefineNodes()); - for (i = 0; i < model.getFEModel()[0]->getNoRefineNodes(); i++) + errors.reserve(nNodes); + for (i = 0; i < nNodes; i++) errors.push_back(DblIdx(0.0,i)); } else { errors.reserve(model.getNoNodes()); @@ -363,7 +372,8 @@ bool AdaptiveSIM::adaptMesh (int iStep) } for (ASMbase* patch : model.getFEModel()) { - if (!patch) return false; + if (!patch) + return false; // extract element norms for this patch Vector locNorm(patch->getNoElms()); @@ -371,14 +381,17 @@ bool AdaptiveSIM::adaptMesh (int iStep) locNorm(i) = eNorm(eRow, patch->getElmID(i)); // remap from geometry basis to refinement basis - Vector locErr(patch->getNoProjectionNodes()); + Vector locErr(patch->getNoRefineNodes()); static_cast(patch)->remapErrors(locErr, locNorm); // insert into global error array for (i = 0; i < locErr.size(); ++i) - if (model.getNoPatches() > 1) - errors[patch->getNodeID(i+1)-1].first += locErr[i]; - else + if (model.getNoPatches() > 1) { + if (prm.MLGN.empty()) + errors[patch->getNodeID(i+1)-1].first += locErr[i]; + else + errors[prm.MLGN[patch->idx][i]].first += locErr[i]; + } else errors[i].first += locErr[i]; } } diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C index 50d69736d..115cad353 100644 --- a/src/SIM/SIMbase.C +++ b/src/SIM/SIMbase.C @@ -1697,7 +1697,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(); diff --git a/src/SIM/SIMinput.C b/src/SIM/SIMinput.C index 32132f035..c7c6e0c7d 100644 --- a/src/SIM/SIMinput.C +++ b/src/SIM/SIMinput.C @@ -1125,9 +1125,15 @@ bool SIMinput::refine (const LR::RefineData& prm, { // Extract local indices from the vector of global indices int locId; - for (int k : prm.elements) - if ((locId = myModel[i]->getNodeIndex(k+1)) > 0) - refineIndices[i].insert(locId-1); + for (int k : prm.elements) { + if (!prm.MLGN.empty()) { + auto it = std::find(prm.MLGN[i].begin(), prm.MLGN[i].end(), k); + if (it != prm.MLGN[i].end()) + refineIndices[i].insert(it-prm.MLGN[i].begin()); + } else + if ((locId = myModel[i]->getNodeIndex(k+1)) > 0) + refineIndices[i].insert(locId-1); + } // fetch all boundary nodes covered (may need to pass this to other patches) pch = dynamic_cast(myModel[i]); @@ -1149,13 +1155,27 @@ bool SIMinput::refine (const LR::RefineData& prm, // for all boundary nodes, check if these appear on other patches for (int k : bndry_nodes) { - int globId = pch->getNodeID(k+1); + int globId; + if (!prm.MLGN.empty()) + globId = prm.MLGN[i][k]; + else + globId = pch->getNodeID(k+1); + for (size_t j = 0; j < myModel.size(); j++) - if (j != i && (locId = myModel[j]->getNodeIndex(globId)) > 0) - { - conformingIndices[j].insert(locId-1); - 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 { + auto it = std::find(prm.MLGN[j].begin(), prm.MLGN[j].end(), globId); + if (it != prm.MLGN[j].end()) { + conformingIndices[j].insert(it-prm.MLGN[j].begin()); + conformingIndices[i].insert(k); + } + } } } diff --git a/src/SIM/SIMinput.h b/src/SIM/SIMinput.h index 5aea15ea8..f2a4fb86b 100644 --- a/src/SIM/SIMinput.h +++ b/src/SIM/SIMinput.h @@ -232,6 +232,9 @@ class SIMinput : public SIMbase //! \brief Deserialization support (for simulation restart). virtual bool deSerialize(const std::map&); + //! \brief Obtain a const reference to model topology. + const std::vector& getInterfaces() const { return myInterfaces; } + private: //! \brief Sets initial conditions from a file. //! \param fieldHolder The SIM-object to inject the initial conditions into