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

Added: SIMgeneric::evalPoint(). Speedup: LRspline2D evaluation through Bezier extraction

Browse Source
This commit is contained in:
Kjetil Andre Johannessen
2015-11-30 10:43:09 +01:00
committed by Knut Morten Okstad
parent ec482d4ff7
commit e2a3e5ecc3
9 changed files with 189 additions and 31 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
cmake/Modules/FindIFEMDeps.cmake
View File

@@ -113,18 +113,18 @@ ENDIF(IFEM_USE_SUPERLU OR IFEM_USE_SUPERLU_MT)
IF(IFEM_USE_LRSPLINES)
FIND_PACKAGE(LRSpline)
IF(LRSpline_LIBRARIES AND LRSpline_INCLUDE_DIRS)
IF(LRSPLINE_VERSION_MINOR LESS 4 AND LRSPLINE_VERSION_MAJOR LESS 1)
IF(LRSPLINE_VERSION_MINOR LESS 5 AND LRSPLINE_VERSION_MAJOR LESS 1)
MESSAGE(STATUS "LR-spline library is out of date.
Found ${LRSPLINE_VERSION_MAJOR}.${LRSPLINE_VERSION_MINOR}.${LRSPLINE_VERSION_PATCH}, need at least 0.4.0.
Found ${LRSPLINE_VERSION_MAJOR}.${LRSPLINE_VERSION_MINOR}.${LRSPLINE_VERSION_PATCH}, need at least 0.5.0.
Support not enabled")
SET(LRSpline_LIBRARIES "")
SET(LRSpline_INCLUDE_DIRS "")
ELSE(LRSPLINE_VERSION_MINOR LESS 4 AND LRSPLINE_VERSION_MAJOR LESS 1)
ELSE(LRSPLINE_VERSION_MINOR LESS 5 AND LRSPLINE_VERSION_MAJOR LESS 1)
SET(IFEM_DEPLIBS ${IFEM_DEPLIBS} ${LRSpline_LIBRARIES})
SET(IFEM_DEPINCLUDES ${IFEM_DEPINCLUDES} ${LRSpline_INCLUDE_DIRS})
SET(IFEM_BUILD_CXX_FLAGS "${IFEM_BUILD_CXX_FLAGS} -DHAS_LRSPLINE=1 ${LRSpline_DEFINITIONS}")
SET(IFEM_CXX_FLAGS "${IFEM_CXX_FLAGS} -DHAS_LRSPLINE=1 ${LRSpline_DEFINITIONS}")
ENDIF(LRSPLINE_VERSION_MINOR LESS 4 AND LRSPLINE_VERSION_MAJOR LESS 1)
ENDIF(LRSPLINE_VERSION_MINOR LESS 5 AND LRSPLINE_VERSION_MAJOR LESS 1)
ENDIF(LRSpline_LIBRARIES AND LRSpline_INCLUDE_DIRS)
ENDIF(IFEM_USE_LRSPLINES)

15
src/ASM/ASMunstruct.C
View File

@@ -26,7 +26,7 @@ int ASMunstruct::gNod = 0;
ASMunstruct::ASMunstruct (unsigned char n_p, unsigned char n_s,
unsigned char n_f)
unsigned char n_f)
: ASMbase(n_p,n_s,n_f)
{
geo = NULL;
@@ -53,15 +53,24 @@ ASMunstruct::~ASMunstruct ()
bool ASMunstruct::refine (const RealArray&, const IntVec&, Vectors*, const char*)
{
std::cerr <<" *** ASMunstruct::refine: Not available without LR-Spline"
<< std::endl;
<< std::endl;
return false;
}
bool ASMunstruct::refine (const IntVec&, const IntVec&, Vectors*, const char*)
{
std::cerr <<" *** ASMunstruct::refine: Not available without LR-Spline"
<< std::endl;
<< std::endl;
return false;
}
Go::BsplineBasis ASMunstruct::getBezierBasis (int p)
{
std::vector<double> knot(2*p);
std::fill(knot.begin(), knot.begin()+p, -1.0);
std::fill(knot.begin()+p, knot.end() , 1.0);
return Go::BsplineBasis(p,p,knot.begin());
}
#endif

6
src/ASM/ASMunstruct.h
View File

@@ -15,6 +15,7 @@
#define _ASM_UNSTRUCT_H
#include "ASMbase.h"
#include "GoTools/geometry/BsplineBasis.h"
namespace LR {
class LRSpline;
@@ -79,6 +80,11 @@ public:
//! \param[in] integrand Object with problem-specific data and methods
virtual LR::LRSpline* evalSolution(const IntegrandBase& integrand) const = 0;
/*!
\brief Returns a Bezier basis of order \a p.
*/
static Go::BsplineBasis getBezierBasis (int p);
protected:
LR::LRSpline* geo; //!< Pointer to the actual spline geometry object

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

@@ -126,8 +126,10 @@ bool ASMunstruct::refine (const RealArray& elementError,
bool isLinIndep = geo->isLinearIndepByOverloading(false);
if(!isLinIndep) {
std::cout << "Inconclusive..." << std::endl;
#ifdef HAS_BOOST
std::cout << "Testing for linear independence by full tensor expansion " << std::endl;
isLinIndep = geo->isLinearIndepByMappingMatrix(false);
#endif
}
if(!isLinIndep)
{
@@ -259,8 +261,10 @@ bool ASMunstruct::refine (const IntVec& elements,
bool isLinIndep = geo->isLinearIndepByOverloading(false);
if(!isLinIndep) {
std::cout << "Inconclusive..." << std::endl;
#ifdef HAS_BOOST
std::cout << "Testing for linear independence by full tensor expansion " << std::endl;
isLinIndep = geo->isLinearIndepByMappingMatrix(false);
#endif
}
if(!isLinIndep)
{
@@ -286,3 +290,11 @@ bool ASMunstruct::refine (const IntVec& elements,
return true;
}
Go::BsplineBasis ASMunstruct::getBezierBasis (int p)
{
double knot[2*p];
std::fill(knot, knot+p, -1.0);
std::fill(knot+p, knot+2*p, 1.0);
return Go::BsplineBasis(p,p,knot);
}

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

@@ -297,6 +297,82 @@ bool ASMu2D::raiseOrder (int ru, int rv)
}
void ASMu2D::evaluateBasis(FiniteElement &fe, int derivs) const
{
// skipping all error testing. Its probably really really smart
LR::Element *el = lrspline->getElement(fe.iel-1);
fe.xi = 2.0*(fe.u - el->umin()) / (el->umax() - el->umin()) - 1.0; // sets value between -1 and 1
fe.eta = 2.0*(fe.v - el->vmin()) / (el->vmax() - el->vmin()) - 1.0;
std::vector<double> Nu = bezier_u.computeBasisValues(fe.xi, derivs);
std::vector<double> Nv = bezier_v.computeBasisValues(fe.eta, derivs);
const Matrix &C = bezierExtract[fe.iel-1];
int p1 = lrspline->order(0);
int p2 = lrspline->order(1);
int p = p1*p2;
Vector B(p); // Bezier basis functions (vector of p1 x p2 components)
Vector Bu(p); // Bezier basis functions differentiated wrt u
Vector Bv(p); // Bezier basis functions differentiated wrt v
size_t k = 0;
for(size_t j=0; j<Nv.size(); j+=(derivs+1)) {
for(size_t i=0; i<Nu.size(); i+=(derivs+1), k++) {
B[k] = Nu[i]*Nv[j];
if(derivs > 0) {
Bu[k] = Nu[i+1]*Nv[ j ];
Bv[k] = Nu[ i ]*Nv[j+1];
}
}
}
fe.N = C*B;
int N = el->nBasisFunctions();
int allP = p1*p2;
double sum = 0;
for(int qq=1; qq<=N; qq++) sum+= fe.N(qq);
if(fabs(sum-1) > 1e-10) {
std::cerr << "fe.N not sums to one at integration point #" << 1 << std::endl;
exit(123);
}
sum = 0;
for(int qq=1; qq<=allP; qq++) sum+= B(qq);
if(fabs(sum-1) > 1e-10) {
std::cerr << "Bezier basis not sums to one at integration point #" << 1 << std::endl;
exit(123);
}
if(derivs > 0) {
Matrix dNdu(el->nBasisFunctions(), 2);
dNdu.fillColumn(1, C*Bu);
dNdu.fillColumn(2, C*Bv);
Matrix Xnod, Jac;
getElementCoordinates(Xnod, fe.iel);
fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu);
sum = 0;
for(int qq=1; qq<=N; qq++) sum+= dNdu(qq,1);
if(fabs(sum) > 1e-10) {
std::cerr << "dNdu not sums to zero at integration point #" << 1 << std::endl;
exit(123);
}
sum = 0;
for(int qq=1; qq<=N; qq++) sum+= dNdu(qq,2);
if(fabs(sum) > 1e-10) {
std::cerr << "dNdv not sums to zero at integration point #" << 1 << std::endl;
exit(123);
}
sum = 0;
for(int qq=1; qq<=allP; qq++) sum+= Bu(qq);
if(fabs(sum) > 1e-10) {
std::cerr << "Bezier derivatives not sums to zero at integration point #" << 1 << std::endl;
exit(123);
}
}
}
bool ASMu2D::generateFEMTopology ()
{
// At this point we are through with the tensor spline object.
@@ -322,24 +398,52 @@ bool ASMu2D::generateFEMTopology ()
else if (shareFE)
return true;
const int p1 = lrspline->order(0);
const int p2 = lrspline->order(1);
myMLGN.resize(nBasis);
myMLGE.resize(nElements);
myMNPC.resize(nElements);
bezierExtract.resize(nElements);
lrspline->generateIDs();
std::vector<double> extrMat ;
std::vector<LR::Element*>::iterator el_it = lrspline->elementBegin();
for (int iel=0; iel<nElements; iel++, el_it++)
{
LR::Element *el = *el_it;
int nSupportFunctions = el->nBasisFunctions();
myMLGE[iel] = ++gEl; // global element number over all patches
myMNPC[iel].resize((*el_it)->nBasisFunctions());
myMNPC[iel].resize(el->nBasisFunctions());
int lnod = 0;
for (LR::Basisfunction *b : (*el_it)->support())
for (LR::Basisfunction *b : el->support())
myMNPC[iel][lnod++] = b->getId();
{
PROFILE("Bezier extraction");
// get bezier extraction matrix
lrspline->getBezierExtraction(iel, extrMat);
int width = p1*p2;
int height = nSupportFunctions;
// wrap the extraction data into a Matrix class
Matrix newM(height, width);
for(int c=1; c<=width; c++)
newM.fillColumn(c, &extrMat[(c-1)*height]);
// keep it for use later
bezierExtract[iel] = newM;
}
}
for (int inod = 0; inod < nBasis; inod++)
myMLGN[inod] = ++gNod;
bezier_u = getBezierBasis(p1);
bezier_v = getBezierBasis(p2);
nnod = gNod;
return true;
@@ -1189,10 +1293,18 @@ int ASMu2D::evalPoint (const double* xi, double* param, Vec3& X) const
param[0] = (1.0-xi[0])*lrspline->startparam(0) + xi[0]*lrspline->endparam(0);
param[1] = (1.0-xi[1])*lrspline->startparam(1) + xi[1]*lrspline->endparam(1);
Go::Point X0;
lrspline->point(X0,param[0],param[1]);
for (unsigned char d = 0; d < nsd; d++)
X[d] = X0[d];
int iel = lrspline->getElementContaining(param[0],param[1]);
FiniteElement fe(lrspline->getElement(iel)->nBasisFunctions());
fe.iel = iel + 1;
fe.u = param[0];
fe.v = param[1];
Matrix Xnod;
getElementCoordinates(Xnod, fe.iel);
evaluateBasis(fe);
X = Xnod * fe.N;
return 0;
}
@@ -1385,6 +1497,10 @@ bool ASMu2D::evalSolution (Matrix& sField, const Vector& locSol,
// Fetch element containing evaluation point.
// Sadly, points are not always ordered in the same way as the elements.
int iel = lrspline->getElementContaining(gpar[0][i],gpar[1][i]);
FiniteElement fe(lrspline->getElement(iel)->nBasisFunctions());
fe.iel = iel + 1;
fe.u = gpar[0][i];
fe.v = gpar[1][i];
utl::gather(MNPC[iel],nComp,locSol,eSol);
// Set up control point (nodal) coordinates for current element
@@ -1396,16 +1512,14 @@ bool ASMu2D::evalSolution (Matrix& sField, const Vector& locSol,
switch (deriv) {
case 0: // Evaluate the solution
lrspline->computeBasis(gpar[0][i],gpar[1][i],spline0,iel);
eSol.multiply(spline0.basisValues,ptSol);
evaluateBasis(fe, deriv);
ptSol = eSol * fe.N;
sField.fillColumn(1+i,ptSol);
break;
case 1: // Evaluate first derivatives of the solution
lrspline->computeBasis(gpar[0][i],gpar[1][i],spline1,iel);
SplineUtils::extractBasis(spline1,ptSol,dNdu);
utl::Jacobian(Jac,dNdX,Xnod,dNdu);
ptDer.multiply(eSol,dNdX);
evaluateBasis(fe, deriv);
ptDer.multiply(eSol,fe.dNdX);
sField.fillColumn(1+i,ptDer);
break;

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

@@ -15,6 +15,7 @@
#define _ASM_U2D_H
#include "ASMunstruct.h"
#include "FiniteElement.h"
#include "ASM2D.h"
#include <memory>
@@ -383,6 +384,9 @@ protected:
//! \param[out] XC Coordinates of the element corners
void getElementCorners(int iel, std::vector<Vec3>& XC) const;
//! \brief Evaluate all basis functions and \a derivs number of derivatives on one element
virtual void evaluateBasis(FiniteElement &el, int derivs=0) const;
public:
//! \brief Returns the number of elements on a boundary.
virtual size_t getNoBoundaryElms(char lIndex, char ldim) const;
@@ -414,6 +418,8 @@ protected:
std::shared_ptr<LR::LRSplineSurface> lrspline; //!< Pointer to the LR-spline surface object
std::vector<Matrix> bezierExtract; //!< Bezier extraction matrices
Go::SplineSurface* tensorspline; //!< Pointer to original tensor spline object
// The tensor spline object is kept for backward compatability with the REFINE
// and RAISEORDER key-words, although we take note that there is a possibility
@@ -422,6 +428,9 @@ protected:
//! Inhomogeneous Dirichlet boundary condition data
std::vector<DirichletEdge> dirich;
Go::BsplineBasis bezier_u; // we're keeping one of these objects to do evaluation on
Go::BsplineBasis bezier_v;
};
#endif

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

@@ -816,18 +816,6 @@ void ASMu3D::getElementCorners (int iEl, Vec3Vec& XC) const
}
}
/*!
\brief Returns a Bezier basis of order \a p.
*/
static Go::BsplineBasis getBezierBasis (int p)
{
double knot[2*p];
std::fill(knot, knot+p, -1.0);
std::fill(knot+p, knot+2*p, 1.0);
return Go::BsplineBasis(p,p,knot);
}
void ASMu3D::evaluateBasis (FiniteElement &el, Matrix &dNdu) const
{

11
src/SIM/SIMgeneric.C
View File

@@ -46,3 +46,14 @@ Vector SIMgeneric::getSolution (const Vector& psol, const double* par,
return tmpVal.getColumn(1);
}
int SIMgeneric::evalPoint(const double* xi, Vec3& X, double* param, int patch) const
{
ASMbase *pch = getPatch(patch-1);
double dummy[2];
if(param == NULL)
return pch->evalPoint(xi, dummy, X);
else
return pch->evalPoint(xi, param, X);
}

9
src/SIM/SIMgeneric.h
View File

@@ -45,6 +45,15 @@ public:
//! \return Evaluated solution values
Vector getSolution(const Vector& psol, const double* par,
int deriv = 0, int patch = 1) const;
//! \brief Evaluates the mapping of the geometry at the given point.
//! \param[in] xi Dimensionless parameters in range [0,1] of the point
//! \param[out] X The Cartesian coordinates of the point
//! \param[out] param The parameters of the point in the knot-span domain
//! \param[in] patch The patch to evaluate
//! \return 0 if the evaluation went good
int evalPoint(const double* xi, Vec3& X,
double* param=NULL, int patch = 1) const;
};
#endif