OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Changed: ASMbase::getBoundaryElms() can optionally return patch-local element indices

Browse Source
This commit is contained in:
Knut Morten Okstad 2023-09-23 23:04:12 +02:00
parent 55184edfef
commit 161a3f3c74
17 changed files with 122 additions and 106 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
src/ASM/ASMbase.C
View File

@ -1630,3 +1630,14 @@ bool ASMbase::writeLagBasis (std::ostream& os, const char* type) const
return true;
}
void ASMbase::getBoundaryElms (int lIndex, IntVec& elms,
int orient, bool local) const
{
this->findBoundaryElms(elms,lIndex,orient);
if (local) return;
for (int& elm : elms)
elm = MLGE[elm]-1;
}

27
src/ASM/ASMbase.h
View File

@ -229,14 +229,15 @@ public:
virtual Vec3 getCoord(size_t inod) const = 0;
//! \brief Returns a matrix with all nodal coordinates within the patch.
//! \param[out] X nsd\f$\times\f$n-matrix, where \a n is the number of nodes
//! \param[in] geo If true returns coordinates for geometry basis
//! in the patch
//! \param[in] geo If \e true, coordinates for the geometry basis are returned
//! otherwise the integration basis coordinates are returned
virtual void getNodalCoordinates(Matrix& X, bool geo = false) const = 0;
//! \brief Returns a matrix with nodal coordinates for an element.
//! \param[in] iel 1-based element index local to current patch
//! \param[in] forceItg If true force returning integration basis coordinates
//! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
//! in one element
//! \param[in] iel 1-based element index local to current patch
//! \param[in] forceItg If \e true, return the integration basis coordinates
//! otherwise the geometry basis coordinates are returned
virtual bool getElementCoordinates(Matrix& X, int iel,
bool forceItg = false) const = 0;
@ -245,7 +246,7 @@ public:
//! \param nodes Array of node numbers
//! \param[in] basis Which basis to grab nodes for (for mixed methods)
//! \param[in] thick Thickness of connection
//! \param[in] orient Local orientation of the boundary face/edge
//! \param[in] orient Local orientation flag (for LR splines only)
//! \param[in] local If \e true, return patch-local numbers
virtual void getBoundaryNodes(int lIndex, IntVec& nodes,
int basis = 0, int thick = 1,
@ -262,11 +263,13 @@ public:
int basis = 0, int orient = -1,
bool local = false, bool open = false) const {}
//! \brief Finds the global (or patch-local) node numbers on a patch boundary.
//! \param[in] lIndex Local index of the boundary face/edge
//! \param[in] orient Local orientation of the boundary face/edge
//! \brief Finds the global (or patch-local) element numbers on a boundary.
//! \param[in] lIndex Local index of the boundary face/edge/vertex
//! \param[out] elms Array of element numbers
virtual void getBoundaryElms(int lIndex, int orient, IntVec& elms) const = 0;
//! \param[in] orient Local orientation flag (for LR splines only)
//! \param[in] local If \e true, return patch-local element numbers
void getBoundaryElms(int lIndex, IntVec& elms,
int orient = -1, bool local = false) const;
//! \brief Returns (1-based) index of a predefined node set in the patch.
virtual int getNodeSetIdx(const std::string&) const { return 0; }
@ -821,6 +824,12 @@ protected:
int searchCtrlPt(RealArray::const_iterator cit, RealArray::const_iterator end,
const Vec3& X, int dimension, double tol = 0.001) const;
//! \brief Finds the patch-local element numbers on a patch boundary.
//! \param[out] elms Array of element numbers
//! \param[in] lIndex Local index of the boundary face/edge/vertex
//! \param[in] orient Local orientation flag (for LR splines only)
virtual void findBoundaryElms(IntVec& elms, int lIndex, int orient) const = 0;
//! \brief Assembles L2-projection matrices for the secondary solution.
//! \param[out] A Left-hand-side matrix
//! \param[out] B Right-hand-side vectors

7
src/ASM/ASMs1D.C
View File

@ -1893,10 +1893,7 @@ void ASMs1D::getElmConnectivities (IntMat& neigh) const
}
void ASMs1D::getBoundaryElms (int lIndex, int, IntVec& elms) const
void ASMs1D::findBoundaryElms (IntVec& elms, int lIndex, int) const
{
if (lIndex == 1)
elms = {0};
else
elms = {MLGE[nel-1]-1};
elms = { lIndex == 1 ? 0 : static_cast<int>(nel)-1 };
}

10
src/ASM/ASMs1D.h
View File

@ -125,11 +125,6 @@ public:
virtual void getBoundaryNodes(int lIndex, IntVec& nodes,
int, int thick, int, bool local) const;
//! \brief Finds the global (or patch-local) node number on a patch end.
//! \param[in] lIndex Local index of the end point
//! \param[out] elms Array of global element numbers
virtual void getBoundaryElms(int lIndex, int, IntVec& elms) const;
//! \brief Finds the node that is closest to the given point.
//! \param[in] X Global coordinates of point to search for
//! \return 1-based nodal index and distance to to the found node
@ -344,6 +339,11 @@ protected:
// Internal utility methods
// ========================
//! \brief Finds the path-local element numbers on a patch boundary.
//! \param[out] elms Array of element numbers
//! \param[in] lIndex Local index of the end point
virtual void findBoundaryElms(IntVec& elms, int lIndex, int = 0) const;
//! \brief Assembles L2-projection matrices for the secondary solution.
//! \param[out] A Left-hand-side matrix
//! \param[out] B Right-hand-side vectors

6
src/ASM/ASMs2D.C
View File

@ -3239,7 +3239,7 @@ void ASMs2D::getElmConnectivities (IntMat& neigh) const
}
void ASMs2D::getBoundaryElms (int lIndex, int, IntVec& elms) const
void ASMs2D::findBoundaryElms (IntVec& elms, int lIndex, int) const
{
int N1m = surf->numCoefs_u() - surf->order_u() + 1;
int N2m = surf->numCoefs_v() - surf->order_v() + 1;
@ -3250,13 +3250,13 @@ void ASMs2D::getBoundaryElms (int lIndex, int, IntVec& elms) const
case 2:
elms.reserve(N2m);
for (int i = 0; i < N2m; ++i)
elms.push_back(MLGE[i*N1m + (lIndex-1)*(N1m-1)] - 1);
elms.push_back(i*N1m + (lIndex-1)*(N1m-1));
break;
case 3:
case 4:
elms.reserve(N1m);
for (int i = 0; i < N1m; ++i)
elms.push_back(MLGE[i + (lIndex-3)*N1m*(N2m-1)] - 1);
elms.push_back(i + (lIndex-3)*N1m*(N2m-1));
}
}

32
src/ASM/ASMs2D.h
View File

@ -80,8 +80,7 @@ protected:
virtual ~BasisFunctionCache() = default;
//! \brief Returns number of elements in each direction.
const std::array<size_t,2>& noElms() const
{ return nel; }
const std::array<size_t,2>& noElms() const { return nel; }
protected:
//! \brief Implementation specific initialization.
@ -105,22 +104,21 @@ protected:
//! \param reduced True to return index for reduced quadrature
size_t index(size_t el, size_t gp, bool reduced) const override;
protected:
//! \brief Configure quadratures.
bool setupQuadrature();
//! \brief Setup integration point parameters.
virtual void setupParameters();
const ASMs2D& patch; //!< Reference to patch cache is for
int basis; //!< Basis to use
std::array<size_t,2> nel{0,0}; //!< Number of elements in each direction
std::array<size_t,2> nel; //!< Number of elements in each direction
private:
//! \brief Obtain structured element indices.
//! \param el Global element index
std::array<size_t,2> elmIndex(size_t el) const;
//! \brief Configure quadratures.
bool setupQuadrature();
};
public:
@ -238,16 +236,17 @@ public:
virtual int getNodeID(size_t inod, bool noAddedNodes = false) const;
//! \brief Returns a matrix with nodal coordinates for an element.
//! \param[in] iel Element index
//! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
//! in one element
//! \param[in] forceItg If true return integration basis element coordinates
//! \param[in] iel 1-based element index local to current patch
//! \param[in] forceItg If \e true, return the integration basis coordinates
//! otherwise the geometry basis coordinates are returned
virtual bool getElementCoordinates(Matrix& X, int iel, bool forceItg = false) const;
//! \brief Returns a matrix with all nodal coordinates within the patch.
//! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
//! \param[in] geo If true return coordinates of geometry basis
//! in the patch
//! \param[in] geo If \e true, coordinates for the geometry basis are returned
//! otherwise the integration basis coordinates are returned
virtual void getNodalCoordinates(Matrix& X, bool geo = false) const;
//! \brief Returns the global coordinates for the given node.
@ -268,11 +267,6 @@ public:
int basis, int thick = 1,
int = 0, bool local = false) const;
//! \brief Finds the global (or patch-local) node numbers on a patch boundary.
//! \param[in] lIndex Local index of the boundary face/edge
//! \param[out] elms Array of element numbers
virtual void getBoundaryElms(int lIndex, int, IntVec& elms) const;
//! \brief Returns the node index for a given corner.
//! \param[in] I -1 or +1 for either umin or umax corner
//! \param[in] J -1 or +1 for either vmin or vmax corner
@ -317,6 +311,7 @@ public:
//! \brief Checks if a separate projection basis is used for this patch.
virtual bool separateProjectionBasis() const;
// Various methods for preprocessing of boundary conditions and patch topology
// ===========================================================================
@ -589,6 +584,11 @@ protected:
// Internal utility methods
// ========================
//! \brief Finds the patch-local element numbers on a patch boundary.
//! \param[out] elms Array of element numbers
//! \param[in] lIndex Local index of the boundary edge
virtual void findBoundaryElms(IntVec& elms, int lIndex, int = 0) const;
//! \brief Assembles L2-projection matrices for the secondary solution.
//! \param[out] A Left-hand-side matrix
//! \param[out] B Right-hand-side vectors

8
src/ASM/ASMs3D.C
View File

@ -3776,7 +3776,7 @@ void ASMs3D::getElmConnectivities (IntMat& neigh) const
}
void ASMs3D::getBoundaryElms (int lIndex, int, IntVec& elms) const
void ASMs3D::findBoundaryElms (IntVec& elms, int lIndex, int) const
{
int N1m = svol->numCoefs(0) - svol->order(0) + 1;
int N2m = svol->numCoefs(1) - svol->order(1) + 1;
@ -3789,21 +3789,21 @@ void ASMs3D::getBoundaryElms (int lIndex, int, IntVec& elms) const
elms.reserve(N2m*N3m);
for (int k = 0; k < N3m; ++k)
for (int j = 0; j < N2m; ++j)
elms.push_back(MLGE[j*N1m + k*N1m*N2m + (lIndex-1)*(N1m-1)] - 1);
elms.push_back(j*N1m + k*N1m*N2m + (lIndex-1)*(N1m-1));
break;
case 3:
case 4:
elms.reserve(N1m*N3m);
for (int k = 0; k < N3m; ++k)
for (int i = 0; i < N1m; ++i)
elms.push_back(MLGE[i + k*N1m*N2m + (lIndex-3)*(N1m*(N2m-1))] - 1);
elms.push_back(i + k*N1m*N2m + (lIndex-3)*(N1m*(N2m-1)));
break;
case 5:
case 6:
elms.reserve(N1m*N2m);
for (int j = 0; j < N2m; ++j)
for (int i = 0; i < N1m; ++i)
elms.push_back(MLGE[i + j*N1m + (lIndex-5)*(N1m*N2m*(N3m-1))] - 1);
elms.push_back(i + j*N1m + (lIndex-5)*(N1m*N2m*(N3m-1)));
}
}

32
src/ASM/ASMs3D.h
View File

@ -96,10 +96,9 @@ protected:
virtual ~BasisFunctionCache() = default;
//! \brief Returns number of elements in each direction.
const std::array<size_t,3>& noElms() const
{ return nel; }
const std::array<size_t,3>& noElms() const { return nel; }
protected:
protected:
//! \brief Implementation specific initialization.
bool internalInit() override;
@ -108,7 +107,7 @@ protected:
//! \brief Calculates basis function info in a single integration point.
//! \param el Element of integration point (0-indexed)
//! \param gp Integratin point on element (0-indexed)
//! \param gp Integration point on element (0-indexed)
//! \param reduced If true, returns values for reduced integration scheme
BasisFunctionVals calculatePt(size_t el, size_t gp, bool reduced) const override;
@ -129,7 +128,7 @@ protected:
int basis; //!< Basis to use
std::array<size_t,3> nel; //!< Number of elements in each direction
private:
private:
//! \brief Obtain structured element indices.
//! \param el Global element index
std::array<size_t,3> elmIndex(size_t el) const;
@ -210,9 +209,9 @@ public:
//! \param[in] dir Parameter direction defining which boundary to return
virtual Go::SplineSurface* getBoundary(int dir, int = 1);
//! \brief Returns the spline volume representing a basis of this patch.
virtual const Go::SplineVolume* getBasis(int basis = 1) const;
//! \brief Returns the spline volume representing a basis of this patch.
virtual Go::SplineVolume* getBasis(int basis = 1);
//! \brief Returns the spline volume representing a basis of this patch.
virtual const Go::SplineVolume* getBasis(int basis = 1) const;
//! \brief Copies the parameter domain from the \a other patch.
virtual void copyParameterDomain(const ASMbase* other);
@ -253,16 +252,17 @@ public:
virtual int getNodeID(size_t inod, bool noAddedNodes = false) const;
//! \brief Returns a matrix with nodal coordinates for an element.
//! \param[in] iel Element index
//! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
//! in one element
//! \param[in] forceItg If true return integration basis element coordinates
//! \param[in] iel 1-based element index local to current patch
//! \param[in] forceItg If \e true, return the integration basis coordinates
//! otherwise the geometry basis coordinates are returned
virtual bool getElementCoordinates(Matrix& X, int iel, bool forceItg = false) const;
//! \brief Returns a matrix with all nodal coordinates within the patch.
//! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
//! \param[in] geo If true returns coordinates for geometry basis
//! in the patch
//! \param[in] geo If \e true, coordinates for the geometry basis are returned
//! otherwise the integration basis coordinates are returned
virtual void getNodalCoordinates(Matrix& X, bool geo = false) const;
//! \brief Returns the global coordinates for the given node.
@ -291,11 +291,6 @@ public:
virtual void getBoundary1Nodes(int lEdge, IntVec& nodes, int basis, int = 0,
bool local = false, bool open = false) const;
//! \brief Finds the global (or patch-local) node numbers on a patch boundary.
//! \param[in] lIndex Local index of the boundary face/edge
//! \param[out] elms Array of element numbers
virtual void getBoundaryElms(int lIndex, int, IntVec& elms) const;
//! \brief Returns the node index for a given corner.
//! \param[in] I -1 or +1 for either umin or umax corner
//! \param[in] J -1 or +1 for either vmin or vmax corner
@ -660,6 +655,11 @@ protected:
// Internal utility methods
// ========================
//! \brief Finds the patch-local element numbers on a patch boundary.
//! \param[out] elms Array of element numbers
//! \param[in] lIndex Local index of the boundary face
virtual void findBoundaryElms(IntVec& elms, int lIndex, int = 0) const;
//! \brief Assembles L2-projection matrices for the secondary solution.
//! \param[out] A Left-hand-side matrix
//! \param[out] B Right-hand-side vectors

2
src/ASM/ASMsupel.h
View File

@ -76,7 +76,7 @@ public:
int, int, int, bool local) const;
//! \brief Dummy method doing nothing.
virtual void getBoundaryElms(int, int, IntVec&) const {}
virtual void findBoundaryElms(IntVec&, int, int) const {}
//! \brief Dummy method doing nothing.
virtual bool getParameterDomain(Real2DMat&, IntVec*) const { return false; }
//! \brief Dummy method doing nothing.

22
src/ASM/LR/ASMu2D.C
View File

@ -2928,7 +2928,7 @@ void ASMu2D::getElmConnectivities (IntMat& neigh) const
}
void ASMu2D::getBoundaryElms (int lIndex, int orient, IntVec& elms) const
void ASMu2D::findBoundaryElms (IntVec& elms, int lIndex, int orient) const
{
std::vector<LR::Element*> elements;
switch (lIndex) {
@ -2939,18 +2939,18 @@ void ASMu2D::getBoundaryElms (int lIndex, int orient, IntVec& elms) const
default: return;
}
// Lambda function for sorting wrt. element centre coordinate
auto&& onMidPoint = [orient,lIndex](LR::Element* a, LR::Element* b)
{
int index = lIndex < 3 ? 1 : 0;
double am = a->midpoint()[index];
double bm = b->midpoint()[index];
return orient == 1 ? bm < am : am < bm;
};
if (orient >= 0)
std::sort(elements.begin(), elements.end(),
[orient,lIndex](LR::Element* a, LR::Element* b)
{
int index = lIndex < 3 ? 1 : 0;
double am = a->midpoint()[index];
double bm = b->midpoint()[index];
return orient == 1 ? bm < am : am < bm;
});
std::sort(elements.begin(),elements.end(),onMidPoint);
for (const LR::Element* elem : elements)
elms.push_back(MLGE[elem->getId()]-1);
elms.push_back(elem->getId());
}

12
src/ASM/LR/ASMu2D.h
View File

@ -190,12 +190,6 @@ public:
virtual void getBoundaryNodes(int lIndex, IntVec& nodes, int basis, int = 1,
int orient = 0, bool local = false) const;
//! \brief Finds the global (or patch-local) node numbers on a patch boundary.
//! \param[in] lIndex Local index of the boundary face/edge
//! \param[in] orient Orientation of boundary (used for sorting)
//! \param[out] elms Array of element numbers
virtual void getBoundaryElms(int lIndex, int orient, IntVec& elms) const;
//! \brief Returns the polynomial order in each parameter direction.
//! \param[out] p1 Order in first (u) direction
//! \param[out] p2 Order in second (v) direction
@ -560,6 +554,12 @@ protected:
// Internal utility methods
// ========================
//! \brief Finds the patch-local element numbers on a patch boundary.
//! \param[out] elms Array of element numbers
//! \param[in] lIndex Local index of the boundary edge
//! \param[in] orient Orientation of boundary (used for sorting)
virtual void findBoundaryElms(IntVec& elms, int lIndex, int orient) const;
//! \brief Assembles L2-projection matrices for the secondary solution.
//! \param[out] A Left-hand-side matrix
//! \param[out] B Right-hand-side vectors

37
src/ASM/LR/ASMu3D.C
View File

@ -2355,33 +2355,32 @@ void ASMu3D::getElmConnectivities (IntMat& neigh) const
}
void ASMu3D::getBoundaryElms (int lIndex, int orient, IntVec& elms) const
void ASMu3D::findBoundaryElms (IntVec& elms, int lIndex, int orient) const
{
std::vector<LR::Element*> elements;
this->getBasis(1)->getEdgeElements(elements,getFaceEnum(lIndex));
std::sort(elements.begin(), elements.end(),
[lIndex,orient](const LR::Element* a, const LR::Element* b)
{
int dir = (lIndex - 1) / 2;
int u = dir == 0 ? 1 : 0;
int v = 1 + (dir != 2 ? 1 : 0);
int idx = (orient & 4) ? v : u;
auto A = a->midpoint();
auto B = b->midpoint();
if (A[idx] != B[idx])
return (orient & 2) ? A[idx] > B[idx] : A[idx] < B[idx];
if (orient >= 0)
std::sort(elements.begin(), elements.end(),
[lIndex,orient](const LR::Element* a, const LR::Element* b)
{
int u = lIndex <= 2 ? 1 : 0;
int v = lIndex >= 5 ? 1 : 2;
int idx = (orient & 4) ? v : u;
std::vector<double> A = a->midpoint();
std::vector<double> B = b->midpoint();
if (A[idx] != B[idx])
return (orient & 2) ? A[idx] > B[idx] : A[idx] < B[idx];
idx = (orient & 4) ? u : v;
if (A[idx] != B[idx])
return (orient & 1) ? A[idx] > B[idx] : A[idx] < B[idx];
return false;
});
idx = (orient & 4) ? u : v;
if (A[idx] != B[idx])
return (orient & 1) ? A[idx] > B[idx] : A[idx] < B[idx];
return false;
});
for (const LR::Element* elem : elements)
elms.push_back(MLGE[elem->getId()]-1);
elms.push_back(elem->getId());
}

12
src/ASM/LR/ASMu3D.h
View File

@ -189,12 +189,6 @@ public:
int orient, bool local,
bool open = false) const;
//! \brief Finds the global (or patch-local) node numbers on a patch boundary.
//! \param[in] lIndex Local index of the boundary face/edge
//! \param[in] orient Orientation of boundary (used for sorting)
//! \param[out] elms Array of element numbers
virtual void getBoundaryElms(int lIndex, int orient, IntVec& elms) const;
//! \brief Returns the node index for a given corner.
//! \param[in] I -1 or +1 for either umin or umax corner
//! \param[in] J -1 or +1 for either vmin or vmax corner
@ -572,6 +566,12 @@ protected:
// Internal utility methods
// ========================
//! \brief Finds the patch-local element numbers on a patch boundary.
//! \param[out] elms Array of element numbers
//! \param[in] lIndex Local index of the boundary face
//! \param[in] orient Orientation of boundary (used for sorting)
virtual void findBoundaryElms(IntVec& elms, int lIndex, int orient) const;
//! \brief Assembles L2-projection matrices for the secondary solution.
//! \param[out] A Left-hand-side matrix
//! \param[out] B Right-hand-side vectors

2
src/ASM/Test/TestASMs1D.C
View File

@ -134,7 +134,7 @@ TEST(TestASMs1D, BoundaryElements)
IntMat neigh(3);
std::array<IntVec,2> n;
for (size_t i = 1; i <= 2; ++i) {
pch1.getBoundaryElms(i, 0, n[i-1]);
pch1.getBoundaryElms(i, n[i-1]);
ASSERT_EQ(n[i-1].size(), 1U);
EXPECT_EQ(n[i-1][0], i == 1 ? 0 : 2);
}

2
src/ASM/Test/TestASMs2D.C
View File

@ -50,7 +50,7 @@ TEST(TestASMs2D, BoundaryElements)
std::array<IntVec,4> n;
for (size_t i = 1; i <= 4; ++i) {
pch1.getBoundaryElms(i, 0, n[i-1]);
pch1.getBoundaryElms(i, n[i-1]);
ASSERT_EQ(n[i-1].size(), ref[i-1].size());
for (size_t j = 0; j < ref[i-1].size(); ++j)
EXPECT_EQ(n[i-1][j], ref[i-1][j]);

2
src/ASM/Test/TestASMs3D.C
View File

@ -64,7 +64,7 @@ TEST(TestASMs3D, BoundaryElements)
std::array<IntVec,6> n;
for (size_t i = 1; i <= 6; ++i) {
pch1.getBoundaryElms(i, 0, n[i-1]);
pch1.getBoundaryElms(i, n[i-1]);
ASSERT_EQ(n[i-1].size(), ref[i-1].size());
for (size_t j = 0; j < ref[i-1].size(); ++j)
EXPECT_EQ(n[i-1][j], ref[i-1][j]);

4
src/SIM/SIMinput.C
View File

@ -1854,8 +1854,8 @@ IntMat SIMinput::getElmConnectivities () const
if (iface.dim == static_cast<int>(nsd)-1)
{
IntVec sElms, mElms;
myModel[iface.slave-1]->getBoundaryElms(iface.sidx, iface.orient, sElms);
myModel[iface.master-1]->getBoundaryElms(iface.midx, 0, mElms);
myModel[iface.slave-1]->getBoundaryElms(iface.sidx, sElms, iface.orient);
myModel[iface.master-1]->getBoundaryElms(iface.midx, mElms, 0);
DomainDecomposition::OrientIterator iter(myModel[iface.slave-1],
iface.orient, iface.sidx);