Skip to content

Commit 4a7099e

Browse files
akva2akva2
committed
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.
1 parent 3b83047 commit 4a7099e

File tree

11 files changed

+307
-31
lines changed

11 files changed

+307
-31
lines changed

‎src/ASM/ASMunstruct.h‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -46,6 +46,7 @@ namespace LR //! Utilities for LR-splines.
4646
IntVec options; //!< Parameters used to control the refinement
4747
IntVec elements; //!< 0-based indices of the elements to refine
4848
RealArray errors; //!< List of error indicators for the elements
49+
std::vector<IntVec> MLGN; //!< MLGN mapping to use for multipatch
4950

5051
//! \brief Default constructor.
5152
explicit RefineData(bool rs = false) : refShare(rs) {}
@@ -214,6 +215,9 @@ class ASMunstruct : public ASMbase
214215
//! \param[in] refTol Mesh refinement threshold
215216
virtual bool refine(const RealFunc& refC, double refTol) = 0;
216217

218+
//! \brief Obtain the refinement basis.
219+
virtual const LR::LRSpline* getRefinementBasis() const { return geo; }
220+
217221
protected:
218222
//! \brief Refines the mesh adaptively.
219223
//! \param[in] prm Input data used to control the mesh refinement

‎src/ASM/GlobalNodes.C‎

Lines changed: 174 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,174 @@
1+
// $Id$
2+
//==============================================================================
3+
//!
4+
//! \file GlobalNodes.C
5+
//!
6+
//! \date Mar 13 2018
7+
//!
8+
//! \author Arne Morten Kvarving / SINTEF
9+
//!
10+
//! \brief Simple global node establishment for unstructured FE models.
11+
//!
12+
//==============================================================================
13+
14+
#include "GlobalNodes.h"
15+
#include "ASMunstruct.h"
16+
#include "Utilities.h"
17+
18+
19+
GlobalNodes::IntVec GlobalNodes::getBoundaryNodes(const LR::LRSpline& lr,
20+
int dim, int lidx, int orient)
21+
{
22+
GlobalNodes::IntVec lNodes;
23+
std::vector<LR::Basisfunction*> edgeFunctions;
24+
LR::parameterEdge edge;
25+
if (dim == 0) {
26+
if (lr.nVariate() == 2) {
27+
switch (lidx) {
28+
case 1: edge = LR::WEST | LR::SOUTH; break;
29+
case 2: edge = LR::EAST | LR::SOUTH; break;
30+
case 3: edge = LR::WEST | LR::NORTH; break;
31+
case 4: edge = LR::EAST | LR::NORTH; break;
32+
}
33+
} else {
34+
switch (lidx) {
35+
case 1: edge = LR::WEST | LR::SOUTH | LR::BOTTOM; break;
36+
case 2: edge = LR::EAST | LR::SOUTH | LR::BOTTOM; break;
37+
case 3: edge = LR::WEST | LR::NORTH | LR::BOTTOM; break;
38+
case 4: edge = LR::EAST | LR::NORTH | LR::BOTTOM; break;
39+
case 5: edge = LR::WEST | LR::SOUTH | LR::TOP; break;
40+
case 6: edge = LR::EAST | LR::SOUTH | LR::TOP; break;
41+
case 7: edge = LR::WEST | LR::NORTH | LR::TOP; break;
42+
case 8: edge = LR::EAST | LR::NORTH | LR::TOP; break;
43+
}
44+
}
45+
} else if (dim == 1) {
46+
if (lr.nVariate() == 2) {
47+
switch (lidx) {
48+
case 1: edge = LR::WEST; break;
49+
case 2: edge = LR::EAST; break;
50+
case 3: edge = LR::SOUTH; break;
51+
case 4: edge = LR::NORTH; break;
52+
default: break;
53+
}
54+
} else {
55+
switch (lidx) {
56+
case 1: edge = LR::BOTTOM | LR::SOUTH; break;
57+
case 2: edge = LR::BOTTOM | LR::NORTH; break;
58+
case 3: edge = LR::TOP | LR::SOUTH; break;
59+
case 4: edge = LR::TOP | LR::NORTH; break;
60+
case 5: edge = LR::BOTTOM | LR::WEST; break;
61+
case 6: edge = LR::BOTTOM | LR::EAST; break;
62+
case 7: edge = LR::TOP | LR::WEST; break;
63+
case 8: edge = LR::TOP | LR::WEST; break;
64+
case 9: edge = LR::SOUTH | LR::WEST; break;
65+
case 10: edge = LR::SOUTH | LR::EAST; break;
66+
case 11: edge = LR::NORTH | LR::WEST; break;
67+
case 12: edge = LR::NORTH | LR::EAST; break;
68+
}
69+
}
70+
} else if (dim == 2) {
71+
switch (lidx) {
72+
case 1: edge = LR::WEST; break;
73+
case 2: edge = LR::EAST; break;
74+
case 3: edge = LR::SOUTH; break;
75+
case 4: edge = LR::NORTH; break;
76+
case 5: edge = LR::BOTTOM; break;
77+
case 6: edge = LR::TOP; break;
78+
}
79+
}
80+
81+
lr.getEdgeFunctions(edgeFunctions, edge);
82+
83+
if (dim == 1) {
84+
if (lr.nVariate() == 2) {
85+
int v = (lidx == 1 || lidx == 2) ? 0 : 1;
86+
int u = 1-v;
87+
ASMunstruct::Sort(u, v, orient, edgeFunctions);
88+
} else {
89+
int dir = (lidx-1)/4;
90+
int u = dir == 0;
91+
int v = 1 + (dir != 2);
92+
ASMunstruct::Sort(u, v, orient, edgeFunctions);
93+
}
94+
} else if (dim == 2) {
95+
int dir = (lidx-1)/2;
96+
int u = dir == 0;
97+
int v = 1 + (dir != 2);
98+
ASMunstruct::Sort(u, v, orient, edgeFunctions);
99+
}
100+
101+
lNodes.reserve(edgeFunctions.size());
102+
for (const LR::Basisfunction* func : edgeFunctions)
103+
lNodes.push_back(func->getId());
104+
105+
return lNodes;
106+
}
107+
108+
109+
class InterfaceOrder {
110+
public:
111+
bool operator()(const ASM::Interface& A, const ASM::Interface& B) const
112+
{
113+
if (A.master != B.master)
114+
return A.master < B.master;
115+
116+
if (A.slave != B.slave)
117+
return A.slave < B.slave;
118+
119+
if (A.dim != B.dim)
120+
return A.dim < B.dim;
121+
122+
return A.midx < B.midx;
123+
}
124+
};
125+
126+
127+
std::vector<GlobalNodes::IntVec>
128+
GlobalNodes::calcGlobalNodes(const GlobalNodes::LRSplineVec& pchs,
129+
const GlobalNodes::InterfaceVec& interfaces)
130+
{
131+
// count total number of nodes
132+
size_t nNodes = 0;
133+
std::vector<GlobalNodes::IntVec> result(pchs.size());
134+
auto it = result.begin();
135+
for (const LR::LRSpline* pch : pchs) {
136+
it->resize(pch->nBasisFunctions());
137+
std::iota(it->begin(), it->end(), nNodes);
138+
nNodes += pch->nBasisFunctions();
139+
++it;
140+
}
141+
142+
// remap common nodes
143+
InterfaceOrder ifOrder;
144+
std::set<ASM::Interface, InterfaceOrder> ifset(ifOrder);
145+
for (const ASM::Interface& it : interfaces)
146+
ifset.insert(it);
147+
for (size_t i = 0; i < pchs.size(); ++i) {
148+
std::map<int,int> old2new;
149+
for (const ASM::Interface& it : ifset) {
150+
if (it.master != (int)i+1)
151+
continue;
152+
153+
IntVec mNodes = getBoundaryNodes(*pchs[i], it.dim, it.midx, 0);
154+
IntVec sNodes = getBoundaryNodes(*pchs[it.slave-1], it.dim, it.sidx, it.orient);
155+
for (size_t n = 0; n < mNodes.size(); ++n)
156+
old2new[result[it.slave-1][sNodes[n]]] = result[i][mNodes[n]];
157+
}
158+
159+
// renumber
160+
for (size_t j = i; j < pchs.size(); ++j)
161+
for (int& it : result[j])
162+
utl::renumber(it, old2new, false);
163+
164+
// compress
165+
if (i > 0) {
166+
int maxNode = *std::max_element(result[i-1].begin(), result[i-1].end());
167+
for (int& n : result[i])
168+
if (n > maxNode)
169+
n = ++maxNode;
170+
}
171+
}
172+
173+
return result;
174+
}

‎src/ASM/GlobalNodes.h‎

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,50 @@
1+
// $Id$
2+
//==============================================================================
3+
//!
4+
//! \file GlobalNodes.h
5+
//!
6+
//! \date Mar 13 2018
7+
//!
8+
//! \author Arne Morten Kvarving / SINTEF
9+
//!
10+
//! \brief Simple global node establishment for unstructured FE models.
11+
//!
12+
//==============================================================================
13+
14+
#ifndef _GLOBAL_NODES_H_
15+
#define _GLOBAL_NODES_H_
16+
17+
#include "Interface.h"
18+
#include <LRSpline/LRSplineSurface.h>
19+
20+
#include <vector>
21+
22+
23+
/*!
24+
\brief Class establishing global node numbers for unstructed FE models.
25+
*/
26+
27+
class GlobalNodes
28+
{
29+
public:
30+
typedef std::vector<int> IntVec; //!< Convenience typedef
31+
typedef std::vector<const LR::LRSpline*> LRSplineVec; //!< Convenience typedef
32+
typedef std::vector<ASM::Interface> InterfaceVec; //!< Convenience typedef
33+
34+
//! \brief Extract local boundary nodes for a LR spline.
35+
//! \param lr The LR spline to extract boundary nodes for
36+
//! \param dim The dimension of the boundary to extract
37+
//! \param lidx The local index of the boundary to extract
38+
//! \param orient Orientation of nodes on boundary
39+
static IntVec getBoundaryNodes(const LR::LRSpline& lr,
40+
int dim, int lidx, int orient);
41+
42+
43+
//! \brief Calculate global node numbers for a FE model.
44+
//! \param pchs The spline patches in the model
45+
//! \param interfaces The topological connections for the spline patches
46+
static std::vector<IntVec> calcGlobalNodes(const LRSplineVec& pchs,
47+
const InterfaceVec& interfaces);
48+
};
49+
50+
#endif

‎src/ASM/LR/ASMLRSpline.C‎

Lines changed: 9 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@
1515
#include "IFEM.h"
1616
#include "LRSpline/LRSplineSurface.h"
1717
#include "LRSpline/Basisfunction.h"
18+
#include "GlobalNodes.h"
1819
#include "Profiler.h"
1920
#include "ThreadGroups.h"
2021
#include <fstream>
@@ -334,15 +335,16 @@ void ASMunstruct::getFunctionsForElements (IntSet& functions,
334335
IntVec ASMunstruct::getBoundaryNodesCovered (const IntSet& nodes) const
335336
{
336337
IntSet result;
338+
const LR::LRSpline* refBasis = this->getRefinementBasis();
337339
int numbEdges = (this->getNoParamDim() == 2) ? 4 : 6;
338340
for (int edge = 1; edge <= numbEdges; edge++)
339341
{
340-
IntVec oneBoundary;
341-
this->getBoundaryNodes(edge, oneBoundary, 1, 1, 0, true); // this returns a 1-indexed list
342+
IntVec bnd = GlobalNodes::getBoundaryNodes(*refBasis,
343+
this->getNoParamDim()-1, edge, 0);
342344
for (const int i : nodes)
343-
for (const int j : oneBoundary)
344-
if (geo->getBasisfunction(i)->contains(*geo->getBasisfunction(j-1)))
345-
result.insert(j-1);
345+
for (const int j : bnd)
346+
if (refBasis->getBasisfunction(i)->contains(*refBasis->getBasisfunction(j)))
347+
result.insert(j);
346348
}
347349

348350
return IntVec(result.begin(), result.end());
@@ -352,9 +354,10 @@ IntVec ASMunstruct::getBoundaryNodesCovered (const IntSet& nodes) const
352354
IntVec ASMunstruct::getOverlappingNodes (const IntSet& nodes, int dir) const
353355
{
354356
IntSet result;
357+
const LR::LRSpline* refBasis = this->getRefinementBasis();
355358
for (const int i : nodes)
356359
{
357-
LR::Basisfunction *b = geo->getBasisfunction(i);
360+
const LR::Basisfunction *b = refBasis->getBasisfunction(i);
358361
for (auto el : b->support()) // for all elements where *b has support
359362
for (auto basis : el->support()) // for all functions on this element
360363
{

‎src/ASM/LR/ASMu2D.C‎

Lines changed: 8 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@
2222
#include "TimeDomain.h"
2323
#include "FiniteElement.h"
2424
#include "GlobalIntegral.h"
25+
#include "GlobalNodes.h"
2526
#include "LocalIntegral.h"
2627
#include "IntegrandBase.h"
2728
#include "CoordinateMapping.h"
@@ -355,10 +356,12 @@ void ASMu2D::extendRefinementDomain (IntSet& refineIndices,
355356

356357
IntVec bndry0;
357358
std::vector<IntVec> bndry1(4);
358-
for (int i = 0; i < 4; i++)
359+
for (int i = 1; i <= 4; i++)
359360
{
360-
bndry0.push_back(this->getCorner((i%2)*2-1, (i/2)*2-1, 1));
361-
this->getBoundaryNodes(1+i, bndry1[i], 1, 1, 0, true);
361+
bndry0.push_back(GlobalNodes::getBoundaryNodes(*this->getRefinementBasis(),
362+
0, i, 0).front());
363+
bndry1[i-1] = GlobalNodes::getBoundaryNodes(*this->getRefinementBasis(),
364+
1, i, 0);
362365
}
363366

364367
// Add refinement from neighbors
@@ -369,7 +372,7 @@ void ASMu2D::extendRefinementDomain (IntSet& refineIndices,
369372
// Check if node is a corner node,
370373
// compute large extended domain (all directions)
371374
for (int edgeNode : bndry0)
372-
if (edgeNode-1 == j)
375+
if (edgeNode == j)
373376
{
374377
IntVec secondary = this->getOverlappingNodes(j);
375378
refineIndices.insert(secondary.begin(),secondary.end());
@@ -381,7 +384,7 @@ void ASMu2D::extendRefinementDomain (IntSet& refineIndices,
381384
// compute small extended domain (one direction)
382385
for (int edge = 0; edge < 4 && !done_with_this_node; edge++)
383386
for (int edgeNode : bndry1[edge])
384-
if (edgeNode-1 == j)
387+
if (edgeNode == j)
385388
{
386389
IntVec secondary = this->getOverlappingNodes(j);
387390
refineIndices.insert(secondary.begin(),secondary.end());

‎src/ASM/LR/ASMu2Dmx.C‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1200,3 +1200,9 @@ bool ASMu2Dmx::evalProjSolution (Matrix& sField, const Vector& locSol,
12001200

12011201
return true;
12021202
}
1203+
1204+
1205+
const LR::LRSpline* ASMu2Dmx::getRefinementBasis() const
1206+
{
1207+
return refBasis.get();
1208+
}

‎src/ASM/LR/ASMu2Dmx.h‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -222,6 +222,9 @@ class ASMu2Dmx : public ASMu2D, private ASMmxBase
222222
virtual void remapErrors(RealArray& errors,
223223
const RealArray& origErr, bool elemErrors) const;
224224

225+
//! \brief Obtain the refinement basis.
226+
virtual const LR::LRSpline* getRefinementBasis() const;
227+
225228
protected:
226229
//! \brief Assembles L2-projection matrices for the secondary solution.
227230
//! \param[out] A Left-hand-side matrix

0 commit comments

Comments
 (0)