diff --git a/Apps/HDF5toVTx/HDF5toVTx.C b/Apps/HDF5toVTx/HDF5toVTx.C index d0b99f003..93a9eaece 100644 --- a/Apps/HDF5toVTx/HDF5toVTx.C +++ b/Apps/HDF5toVTx/HDF5toVTx.C @@ -80,6 +80,7 @@ bool readBasis (std::vector& result, const std::string& name, std::string out; hdf.readString(geom.str(),out); ptype = out.substr(0,10) == "# LRSPLINE" ? ASM::LRSpline : ASM::Spline; + isLR = ptype == ASM::LRSpline; int gdim = 0; if (ptype == ASM::LRSpline) gdim = out.substr(11,7) == "SURFACE" ? 2 : 3; 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 34a2ce353..df2ba1ae7 100644 --- a/src/ASM/ASMunstruct.h +++ b/src/ASM/ASMunstruct.h @@ -46,6 +46,8 @@ 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 + IntVec pMLGN; //!< Parallel MLGN mapping to use with MPI //! \brief Default constructor. explicit RefineData(bool rs = false) : refShare(rs) {} @@ -217,6 +219,9 @@ class ASMunstruct : public ASMbase //! \brief Finds the node that is closest to the given point \b X. virtual std::pair findClosestNode(const Vec3& X) const; + //! \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/DomainDecomposition.C b/src/ASM/DomainDecomposition.C index bfc344100..bab07a930 100644 --- a/src/ASM/DomainDecomposition.C +++ b/src/ASM/DomainDecomposition.C @@ -322,8 +322,10 @@ void DomainDecomposition::setupNodeNumbers(int basis, IntVec& lNodes, // need to expand to all bases for corners and edges for (size_t b = 1; b <= pch->getNoBasis(); ++b) cbasis.insert(b); - else // directly add nodes, cbasis remains empty + else { // directly add nodes, cbasis remains empty + std::cout << "we go here yo" << std::endl; pch->getBoundaryNodes(lidx, lNodes, 0, thick, orient, false); + } const ASM2D* pch2D = dynamic_cast(pch); const ASM3D* pch3D = dynamic_cast(pch); diff --git a/src/ASM/GlobalNodes.C b/src/ASM/GlobalNodes.C new file mode 100644 index 000000000..aaf11a851 --- /dev/null +++ b/src/ASM/GlobalNodes.C @@ -0,0 +1,259 @@ +// $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 "SIMbase.h" +#include "DomainDecomposition.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; +} + + +/*! + \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; + } +}; + + +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; +} + + +GlobalNodes::IntVec +GlobalNodes::calcDDMapping(const GlobalNodes::LRSplineVec& pchs, + const std::vector& MLGN, + const SIMbase& sim, int& nNodes) +{ + const ProcessAdm& adm = sim.getProcessAdm(); + + int minNode = 0; + if (adm.getProcId() > 0) + adm.receive(minNode, adm.getProcId()-1); + + IntVec result(*std::max_element(MLGN.back().begin(), MLGN.back().end()) + 1); + std::iota(result.begin(), result.end(), minNode); + int maxNode = adm.getProcId() == 0 ? 0 : minNode; + + std::map old2new; + for (const auto& it : adm.dd.ghostConnections) { + int sidx = sim.getLocalPatchIndex(it.slave); + if (sidx < 1) + continue; + + IntVec lNodes = GlobalNodes::getBoundaryNodes(*pchs[sidx-1], it.dim, it.sidx, 0); + for (int& i : lNodes) + i = result[MLGN[sidx-1][i]]; + + int nRecv; + adm.receive(nRecv, adm.dd.getPatchOwner(it.master)); + if (nRecv =! lNodes.size()) { + std::cerr <<"\n *** GlobalNodes::calcDDMapping(): " + <<" Topology error, boundary size " + << nRecv << ", expected " << lNodes.size() << std::endl; + return IntVec(); + } + IntVec glbNodes(lNodes.size()); + adm.receive(glbNodes, adm.dd.getPatchOwner(it.master)); + for (size_t i = 0; i < lNodes.size(); ++i) + old2new[lNodes[i]] = glbNodes[i]; + } + + // remap ghost nodes + for (auto& it : result) + utl::renumber(it, old2new, false); + + // remap rest of our nodes + for (int i = 0; i < (int)result.size(); ++i) + if (old2new.find(i + minNode) == old2new.end()) { + std::map old2new2; + old2new2[i + minNode] = maxNode++; + for (auto& it : result) + utl::renumber(it, old2new2, false); + } + + if (adm.getProcId() < adm.getNoProcs()-1) + adm.send(maxNode, adm.getProcId()+1); + + for (const auto& it : adm.dd.ghostConnections) { + int midx = sim.getLocalPatchIndex(it.master); + if (midx < 1) + continue; + + IntVec glbNodes = GlobalNodes::getBoundaryNodes(*pchs[midx-1], it.dim, + it.midx, it.orient); + for (size_t i = 0; i < glbNodes.size(); ++i) + glbNodes[i] = result[MLGN[midx-1][glbNodes[i]]]; + + adm.send(int(glbNodes.size()), adm.dd.getPatchOwner(it.slave)); + adm.send(glbNodes, adm.dd.getPatchOwner(it.slave)); + } + +#ifdef HAVE_MPI + nNodes = adm.allReduce(adm.getProcId() == adm.getNoProcs()-1 ? maxNode : 0, MPI_SUM); +#endif + + return result; +} diff --git a/src/ASM/GlobalNodes.h b/src/ASM/GlobalNodes.h new file mode 100644 index 000000000..7f95af18e --- /dev/null +++ b/src/ASM/GlobalNodes.h @@ -0,0 +1,64 @@ +// $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 + +class DomainDecomposition; +class ProcessAdm; +class SIMbase; + + +/*! + \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); + + //! \brief Calculate parallel global node numbers for a FE model. + //! \param pchs The spline patches in the model + //! \param MLGN Process-global node numbers + //! \param sim The simulator holding patch owner information + //! \param adm The parallel process administrator + //! \param dd The domain decomposition holding ghost connections + //! \param[out] nNodes Total number of nodes in the model + static IntVec calcDDMapping(const LRSplineVec& pchs, + const std::vector& MLGN, + const SIMbase& sim, int& nNodes); +}; + +#endif diff --git a/src/ASM/LR/ASMLRSpline.C b/src/ASM/LR/ASMLRSpline.C index 4651d0b22..70ad122f0 100644 --- a/src/ASM/LR/ASMLRSpline.C +++ b/src/ASM/LR/ASMLRSpline.C @@ -17,6 +17,7 @@ #include "LRSpline/Basisfunction.h" #include "Vec3.h" #include "Vec3Oper.h" +#include "GlobalNodes.h" #include "ThreadGroups.h" #include "Profiler.h" #include @@ -332,15 +333,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()); @@ -350,9 +352,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 01f7c252f..88eb19053 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()); @@ -1991,10 +1994,17 @@ void ASMu2D::getBoundaryNodes (int lIndex, IntVec& nodes, int basis, default: return; } - if (nodes.empty()) + if (nodes.empty()) { nodes = this->getEdgeNodes(edge, basis, 0); - else { + if (!local) + for (int& node : nodes) + node = MLGN[node-1]; + } else { IntVec nodes2 = this->getEdgeNodes(edge, basis, 0); + if (!local) + for (int& node : nodes2) + node = MLGN[node-1]; + nodes.insert(nodes.end(), nodes2.begin(), nodes2.end()); } diff --git a/src/ASM/LR/ASMu2Dmx.C b/src/ASM/LR/ASMu2Dmx.C index d7e86276f..e71f61a1c 100644 --- a/src/ASM/LR/ASMu2Dmx.C +++ b/src/ASM/LR/ASMu2Dmx.C @@ -1140,10 +1140,7 @@ size_t ASMu2Dmx::getNoProjectionNodes() const Fields* ASMu2Dmx::getProjectedFields(const Vector& coefs, size_t nf) const { - if (projBasis != m_basis[0]) - return new LRSplineFields2D(projBasis.get(), coefs, nf); - - return nullptr; + return new LRSplineFields2D(projBasis.get(), coefs, nf); } @@ -1200,3 +1197,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..9167e7377 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; } @@ -2086,18 +2082,27 @@ void ASMu3D::getBoundaryNodes (int lIndex, IntVec& nodes, int basis, if (!this->getBasis(basis)) return; // silently ignore empty patches - nodes = this->getFaceNodes(lIndex, basis, orient); + if (nodes.empty()) { + std::cout << "herrrr " << basis << std::endl; + nodes = this->getFaceNodes(lIndex, basis, orient); + if (!local) + for (int& node : nodes) + node = this->getNodeID(node); + } else { + IntVec nodes2 = this->getFaceNodes(lIndex, basis, orient); + if (!local) + for (int& node : nodes2) + node = this->getNodeID(node); + nodes.insert(nodes.end(), nodes2.begin(), nodes2.end()); + } -#if SP_DEBUG > 1 +//#if SP_DEBUG > 1 std::cout <<"Boundary nodes in patch "<< idx+1 <<" edge "<< lIndex <<":"; for (size_t i = 0; i < nodes.size(); i++) std::cout <<" "<< nodes[i]; std::cout << std::endl; -#endif +//#endif - if (!local) - for (int& node : nodes) - node = this->getNodeID(node); } @@ -2478,18 +2483,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 +2505,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 +2518,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 +2536,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..703420f09 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]); @@ -1037,10 +1040,7 @@ size_t ASMu3Dmx::getNoRefineElms() const Fields* ASMu3Dmx::getProjectedFields(const Vector& coefs, size_t nf) const { - if (projBasis != m_basis[0]) - return new LRSplineFields3D(projBasis.get(), coefs, nf); - - return nullptr; + return new LRSplineFields3D(projBasis.get(), coefs, nf); } @@ -1086,3 +1086,36 @@ 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(); +} + + +void ASMu3Dmx::getBoundaryNodes (int lIndex, IntVec& nodes, int basis, + int thick, int orient, bool local) const +{ + if (basis > 0) + this->ASMu3D::getBoundaryNodes(lIndex, nodes, basis, thick, orient, local); + else + for (size_t b = 1; b <= this->getNoBasis(); ++b) + this->ASMu3D::getBoundaryNodes(lIndex, nodes, b, thick, orient, local); +} diff --git a/src/ASM/LR/ASMu3Dmx.h b/src/ASM/LR/ASMu3Dmx.h index 3c2345d99..00d833ab9 100644 --- a/src/ASM/LR/ASMu3Dmx.h +++ b/src/ASM/LR/ASMu3Dmx.h @@ -185,6 +185,17 @@ class ASMu3Dmx : public ASMu3D, private ASMmxBase virtual bool refine(const LR::RefineData& prm, Vectors& sol, const char* fName = nullptr); + //! \brief Finds the global (or patch-local) node numbers on a patch boundary. + //! \param[in] lIndex Local index of the boundary edge + //! \param nodes Array of node numbers + //! \param[in] basis Which basis to grab nodes for (0 for all) + //! \param[in] thick Thickness of connection + //! \param[in] orient Orientation for returned boundary nodes + //! \param[in] local If \e true return patch-local node numbers + virtual void getBoundaryNodes(int lIndex, IntVec& nodes, + int basis, int thick = 1, + int orient = 0, bool local = false) const; + //! \brief Remap (geometry) element wise errors to refinement basis functions. //! \param errors The remapped errors //! \param[in] origErr The element wise errors on the geometry mesh @@ -192,6 +203,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 991ecb3f4..41f32d93e 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" @@ -308,9 +311,9 @@ bool AdaptiveSIM::adaptMesh (int iStep, std::streamsize outPrec) if (eNorm.cols() < 1 || gNorm[adaptor](adNorm) < errTol*refNorm) return false; // Calculate row index in eNorm of the error norm to adapt based on - size_t i, eRow = adNorm; + size_t eRow = adNorm; NormBase* norm = model.getNormIntegrand(); - for (i = 0; i < adaptor; i++) + for (size_t i = 0; i < adaptor; i++) eRow += norm->getNoFields(i+1); delete norm; @@ -351,38 +354,47 @@ bool AdaptiveSIM::adaptMesh (int iStep, std::streamsize outPrec) { ASMbase* patch = model.getPatch(1); if (!patch) return false; - - if (patch->getNoRefineNodes() == patch->getNoNodes(1)) - errors.resize(model.getNoNodes(),DblIdx(0.0,0)); - else if (model.getNoPatches() == 1) - errors.resize(patch->getNoRefineNodes(),DblIdx(0.0,0)); - else - { - std::cerr <<" *** AdaptiveSIM::adaptMesh: Multi-patch refinement" - <<" is not available for mixed models."<< std::endl; - return false; + size_t nNodes = patch->getNoRefineNodes(); + if (model.getNoPatches() > 1) { + 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; + if (model.getProcessAdm().getNoProcs() > 1) + prm.pMLGN = GlobalNodes::calcDDMapping(refBasis, prm.MLGN, model, nNodes); +#endif } + errors.reserve(nNodes); + for (int i = 0; i < nNodes; i++) + errors.push_back(DblIdx(0.0,i)); for (i = 0; i < errors.size(); i++) errors[i].second = i; for (ASMbase* patch : model.getFEModel()) { - if (!patch) return false; + if (!patch) + return false; // extract element norms for this patch Vector locNorm(patch->getNoElms()); - for (i = 1; i <= patch->getNoElms(); ++i) + for (size_t i = 1; i <= patch->getNoElms(); ++i) 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 + for (size_t i = 0; i < locErr.size(); ++i) + if (model.getNoPatches() > 1) { + if (prm.pMLGN.empty()) + errors[prm.MLGN[patch->idx][i]].first += locErr[i]; + else + errors[prm.pMLGN[prm.MLGN[patch->idx][i]]].first += locErr[i]; + } else errors[i].first += locErr[i]; } } @@ -393,10 +405,21 @@ bool AdaptiveSIM::adaptMesh (int iStep, std::streamsize outPrec) <<" is available for isotropic_function only."<< std::endl; return false; } - for (i = 0; i < eNorm.cols(); i++) + for (size_t i = 0; i < eNorm.cols(); i++) errors.push_back(DblIdx(eNorm(eRow,1+i),i)); } +#ifdef HAVE_MPI + std::vector perr(errors.size(), 0.0); + if (model.getProcessAdm().getNoProcs() > 1) { + for (size_t i = 0; i < errors.size(); ++i) + perr[i] = errors[i].first; + model.getProcessAdm().allReduce(perr, MPI_SUM); + for (size_t i = 0; i < errors.size(); ++i) + errors[i].first = perr[i]; + } +#endif + // Sort the elements in the sequence of decreasing errors std::sort(errors.begin(),errors.end(),std::greater()); @@ -461,8 +484,15 @@ bool AdaptiveSIM::adaptMesh (int iStep, std::streamsize outPrec) return false; prm.elements.reserve(refineSize); - for (i = 0; i < refineSize; i++) - prm.elements.push_back(errors[i].second); + for (size_t i = 0; i < refineSize; i++) + if (prm.pMLGN.empty()) + prm.elements.push_back(errors[i].second); + else { + auto it = std::find(prm.pMLGN.begin(), + prm.pMLGN.end(), errors[i].second); + if (it != prm.pMLGN.end()) + prm.elements.push_back(it-prm.pMLGN.begin()); + } // Now refine the mesh if (!storeMesh) diff --git a/src/SIM/SIMbase.C b/src/SIM/SIMbase.C index 5897a3583..607ced8aa 100644 --- a/src/SIM/SIMbase.C +++ b/src/SIM/SIMbase.C @@ -99,7 +99,6 @@ void SIMbase::clearProperties () for (auto& i4 : myTracs) delete i4.second; - myPatches.clear(); myGlb2Loc.clear(); myScalars.clear(); myVectors.clear(); @@ -1705,7 +1704,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 22faed3a2..72569c8d4 100644 --- a/src/SIM/SIMinput.C +++ b/src/SIM/SIMinput.C @@ -1124,7 +1124,7 @@ bool SIMinput::refine (const LR::RefineData& prm, } // Single patch models only pass refinement call to the ASM level - if (myModel.size() == 1) + if (this->getNoPatches() == 1) return (isRefined = pch->refine(prm,sol,fName)); if (!prm.errors.empty()) // refinement by true_beta @@ -1142,9 +1142,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]); @@ -1166,15 +1172,72 @@ 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); + } + } + } + } + +#ifdef HAVE_MPI + if (this->adm.getNoProcs() > 1) { + for (int i = 0; i < this->adm.getNoProcs(); ++i) { + int send = 0; + if (i == this->adm.getProcId()) + for (auto& it : conformingIndices) + send += it.size(); + + int nRecv = adm.allReduce(send, MPI_SUM); + std::vector vRecv; + if (i == this->adm.getProcId()) { + std::vector vSend; + vRecv.reserve(nRecv); + for (auto& it : conformingIndices) + for (const int& node : it) + vRecv.push_back(prm.pMLGN[node]); + } else + vRecv.resize(nRecv); + + MPI_Bcast(vRecv.data(), vRecv.size(), MPI_INT, i, *this->adm.getCommunicator()); + + if (i != this->adm.getProcId()) { + for (int& node : vRecv) { + auto it = std::find(prm.pMLGN.begin(), prm.pMLGN.end(), node); + if (it != prm.pMLGN.end()) { + int globId = it - prm.pMLGN.begin(); + int locId; + for (size_t j = 0; j < myModel.size(); j++) + if (prm.MLGN.empty()) { + if ((locId = myModel[j]->getNodeIndex(globId)) > 0) + conformingIndices[j].insert(locId-1); + } 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()); + } + } } + } } } +#endif Vectors lsols; lsols.reserve(sol.size()*myModel.size()); 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