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

Fixed: Override getNoIntPoints and getNoBouPoints

Browse Source 
to reflect the respective integrate methods for LRSplines
This commit is contained in:
Knut Morten Okstad 2021-03-15 08:25:47 +01:00
parent a2cec872cb
commit f3cdb440a2
4 changed files with 92 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

89
src/ASM/LR/ASMLRSpline.C
View File

@ -99,24 +99,23 @@ void LR::generateThreadGroups (ThreadGroups& threadGroups,
#ifdef USE_OPENMP
if (omp_get_max_threads() > 1)
{
for (int i = 0; i < 2; i++)
threadGroups[i].clear();
IntMat& answer = threadGroups[0];
threadGroups[0].clear();
threadGroups[1].clear();
IntVec status(nElement,0); // status vector for elements:
std::vector<std::set<int>> additionals;
std::vector<IntSet> additionals;
if (!addConstraints.empty()) {
additionals.resize(nElement);
for (auto e : lr->getAllElements()) {
for (const LRSpline* lr2 : addConstraints) {
for (LR::Element* e : lr->getAllElements()) {
for (LRSpline* lr2 : addConstraints) {
int elB = lr2->getElementContaining(e->midpoint());
for (auto b2 : lr2->getElement(elB)->support())
for (auto el3 : b2->support()) {
std::vector<double> midpoint = el3->midpoint();
std::vector<std::vector<double>> points;
for (LR::Basisfunction* b2 : lr2->getElement(elB)->support())
for (LR::Element* el3 : b2->support()) {
RealArray midpoint = el3->midpoint();
std::vector<RealArray> points;
if (lr2->nElements() != lr->nElements()) {
std::vector<double> diff(midpoint.size());
RealArray diff(midpoint.size());
for (size_t j = 0; j < midpoint.size(); ++j)
diff[j] = (el3->getParmax(j) - el3->getParmin(j)) / 4.0;
if (lr->dimension() == 2)
@ -136,7 +135,7 @@ void LR::generateThreadGroups (ThreadGroups& threadGroups,
} else
points = {midpoint};
for (const std::vector<double>& vec : points)
for (const RealArray& vec : points)
additionals[e->getId()].insert(lr->getElementContaining(vec));
}
}
@ -156,14 +155,14 @@ void LR::generateThreadGroups (ThreadGroups& threadGroups,
// look for available elements
IntVec thisColor;
for (auto e : lr->getAllElements()) {
for (LR::Element* e : lr->getAllElements()) {
int i = e->getId();
if (status[i] == 0) {
status[i] = nColors+1;
thisColor.push_back(i);
fixedElements++;
for (auto b : e->support()) // for all basisfunctions with support here
for (auto el2 : b->support()) {// for all elements this function supports
for (LR::Basisfunction* b : e->support())
for (LR::Element* el2 : b->support()) {
int j = el2->getId();
if (status[j] == 0) // if not assigned a color yet
status[j] = -1; // set as unavailable (with current color)
@ -174,7 +173,7 @@ void LR::generateThreadGroups (ThreadGroups& threadGroups,
}
}
}
answer.push_back(thisColor);
threadGroups[0].push_back(thisColor);
}
return;
}
@ -472,30 +471,28 @@ Vec3 ASMLRSpline::getElementCenter (int iel) const
bool ASMLRSpline::checkThreadGroups (const IntMat& groups,
const std::vector<const LR::LRSpline*> bases,
const std::vector<const LR::LRSpline*>& bases,
const LR::LRSpline* threadBasis)
{
size_t groupId = 1;
bool ok = true;
for (const IntVec& elms : groups) {
size_t basisId = 1;
for (const LR::LRSpline* basis : bases) {
std::set<int> nodes;
for (int elm : elms) {
std::vector<double> midpoint = threadBasis->getElement(elm)->midpoint();
for (size_t gId = 1; gId <= groups.size(); gId++)
for (size_t bId = 1; bId <= bases.size(); bId++)
{
IntSet nodes;
const LR::LRSpline* basis = bases[bId-1];
for (int elm : groups[gId-1]) {
RealArray midpoint = threadBasis->getElement(elm)->midpoint();
int bElm = basis->getElementContaining(midpoint);
for (const LR::Basisfunction* func : basis->getElement(bElm)->support()) {
if (nodes.find(func->getId()) != nodes.end()) {
std::cerr << " ** Warning: Function " << func->getId() << " on basis " << basisId << " present for multiple elements in group " << groupId << std::endl;
for (const LR::Basisfunction* func : basis->getElement(bElm)->support())
if (!nodes.insert(func->getId()).second) {
std::cerr <<" *** ASMLRSpline::checkThreadGroups: Function "
<< func->getId() <<" on basis "<< bId
<<" is present for multiple elements in group "<< gId
<< std::endl;
ok = false;
}
nodes.insert(func->getId());
}
}
++basisId;
}
++groupId;
}
return ok;
}
@ -559,3 +556,31 @@ bool ASMLRSpline::getParameterDomain (Real2DMat& u, IntVec* corners) const
return true;
}
void ASMLRSpline::getNoIntPoints (size_t& nPt, size_t& nIPt)
{
size_t nGp = 1;
if (nGauss > 0 && nGauss <= 10)
for (unsigned char d = 0; d < ndim; d++)
nGp *= nGauss;
else
{
// Use polynomial order to define number of quadrature points
int ng[3] = { 0, 0, 0 };
int nG1 = 0;
this->getOrder(ng[0],ng[1],ng[2]);
for (unsigned char d = 0; d < ndim; d++)
nG1 = std::max(nG1,ng[d]+nGauss%10);
for (unsigned char d = 0; d < ndim; d++)
nGp = nG1*nG1;
}
firstIp = nPt;
nPt += nel*nGp;
// Count additional interface quadrature points
size_t nInterface = MLGE.size() - nel;
if (nInterface > 0 && nInterface != nel && nGauss > 0 && nGauss <= 10)
nIPt += nInterface*nGp/nGauss;
}

13
src/ASM/LR/ASMLRSpline.h
View File

@ -178,6 +178,9 @@ public:
//! \brief Returns the coordinates of the element center.
virtual Vec3 getElementCenter(int iel) const;
//! \brief Computes the total number of integration points in this patch.
virtual void getNoIntPoints(size_t& nPt, size_t& nIPt);
protected:
//! \brief Refines the mesh adaptively.
//! \param[in] prm Input data used to control the mesh refinement
@ -188,15 +191,15 @@ protected:
virtual int evalPoint(int iel, const double* param, Vec3& X) const = 0;
//! \brief Santity check thread groups.
//! \param groups The generated thread groups
//! \param bases The bases to check for
//! \param threadBasis The LRSpline the element groups are derived from
//! \param[in] groups The generated thread groups
//! \param[in] bases The bases to check for
//! \param[in] threadBasis The LRSpline the element groups are derived from
static bool checkThreadGroups(const IntMat& groups,
const std::vector<const LR::LRSpline*> bases,
const std::vector<const LR::LRSpline*>& bases,
const LR::LRSpline* threadBasis);
//! \brief Analyze and print thread group statistics.
//! \param groups The generated thread groups
//! \param[in] groups The generated thread groups
static void analyzeThreadGroups(const IntMat& groups);
LR::LRSpline* geo; //!< Pointer to the actual spline geometry object

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

@ -1982,6 +1982,30 @@ int ASMu3D::getCorner(int I, int J, int K, int basis) const
}
void ASMu3D::getNoBouPoints (size_t& nPt, char ldim, char lindx)
{
size_t nGp = 1;
if (nGauss > 0 && nGauss <= 10)
for (char d = 0; d < ldim; d++)
nGp *= nGauss;
else if (ldim == 2)
{
// Use polynomial order to define number of quadrature points
int p[3] = { 0, 0, 0 };
this->getOrder(p[0],p[1],p[2]);
p[(lindx-1)/2] = 1;
int nG = std::max(std::max(p[0],p[1]),p[2]);
nGp = nG*nG;
}
else
nGp = 0;
firstBp[lindx] = nPt;
nPt += this->getNoBoundaryElms(lindx,ldim)*nGp;
}
void ASMu3D::generateThreadGroups (const Integrand& integrand, bool silence,
bool ignoreGlobalLM)
{

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

@ -277,6 +277,9 @@ public:
// These are the main computational methods of the ASM class hierarchy.
// ====================================================================
//! \brief Computes the number of boundary integration points in this patch.
virtual void getNoBouPoints(size_t& nPt, char ldim, char lindx);
//! \brief Evaluates an integral over the interior patch domain.
//! \param integrand Object with problem-specific data and methods
//! \param glbInt The integrated quantity