Skip to content

Commit

Permalink
added: class for simple global node number establishment for LR splines
Browse files Browse the repository at this point in the history
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.
  • Loading branch information
akva2 committed Mar 15, 2018
1 parent 3b83047 commit e2f0356
Show file tree
Hide file tree
Showing 12 changed files with 313 additions and 31 deletions.
3 changes: 3 additions & 0 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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})
Expand Down
4 changes: 4 additions & 0 deletions src/ASM/ASMunstruct.h
Original file line number Diff line number Diff line change
Expand Up @@ -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<IntVec> MLGN; //!< MLGN mapping to use for multipatch

//! \brief Default constructor.
explicit RefineData(bool rs = false) : refShare(rs) {}
Expand Down Expand Up @@ -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
Expand Down
174 changes: 174 additions & 0 deletions src/ASM/GlobalNodes.C
Original file line number Diff line number Diff line change
@@ -0,0 +1,174 @@
// $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)
{
GlobalNodes::IntVec lNodes;
std::vector<LR::Basisfunction*> edgeFunctions;
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;
}
}

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);
}

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::IntVec>
GlobalNodes::calcGlobalNodes(const GlobalNodes::LRSplineVec& pchs,
const GlobalNodes::InterfaceVec& interfaces)
{
// count total number of nodes
size_t nNodes = 0;
std::vector<GlobalNodes::IntVec> 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<ASM::Interface, InterfaceOrder> ifset(ifOrder);
for (const ASM::Interface& it : interfaces)
ifset.insert(it);
for (size_t i = 0; i < pchs.size(); ++i) {
std::map<int,int> 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
if (i > 0) {
int maxNode = *std::max_element(result[i-1].begin(), result[i-1].end());
for (int& n : result[i])
if (n > maxNode)
n = ++maxNode;
}
}

return result;
}
50 changes: 50 additions & 0 deletions src/ASM/GlobalNodes.h
Original file line number Diff line number Diff line change
@@ -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 <LRSpline/LRSplineSurface.h>

#include <vector>


/*!
\brief Class establishing global node numbers for unstructed FE models.
*/

class GlobalNodes
{
public:
typedef std::vector<int> IntVec; //!< Convenience typedef
typedef std::vector<const LR::LRSpline*> LRSplineVec; //!< Convenience typedef
typedef std::vector<ASM::Interface> 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<IntVec> calcGlobalNodes(const LRSplineVec& pchs,
const InterfaceVec& interfaces);
};

#endif
15 changes: 9 additions & 6 deletions src/ASM/LR/ASMLRSpline.C
Original file line number Diff line number Diff line change
Expand Up @@ -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 <fstream>
Expand Down Expand Up @@ -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());
Expand All @@ -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
{
Expand Down
13 changes: 8 additions & 5 deletions src/ASM/LR/ASMu2D.C
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down Expand Up @@ -355,10 +356,12 @@ void ASMu2D::extendRefinementDomain (IntSet& refineIndices,

IntVec bndry0;
std::vector<IntVec> 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
Expand All @@ -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());
Expand All @@ -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());
Expand Down
6 changes: 6 additions & 0 deletions src/ASM/LR/ASMu2Dmx.C
Original file line number Diff line number Diff line change
Expand Up @@ -1200,3 +1200,9 @@ bool ASMu2Dmx::evalProjSolution (Matrix& sField, const Vector& locSol,

return true;
}


const LR::LRSpline* ASMu2Dmx::getRefinementBasis() const
{
return refBasis.get();
}
3 changes: 3 additions & 0 deletions src/ASM/LR/ASMu2Dmx.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading

0 comments on commit e2f0356

Please sign in to comment.