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

Added: Define quadrature scheme automatically from polynomial order.

Browse Source 
A non-positive value of nGauss is interpreted such that the number of
Gauss points is equal to p+nGauss, where p is the polynomial order.
Changed: Using FiniteElement::operator<< for debug print.
This commit is contained in:
Knut Morten Okstad 2018-03-22 11:52:28 +01:00
parent 536ab5bfd1
commit 8374c8b293
11 changed files with 264 additions and 185 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
src/ASM/ASMbase.C
View File

@ -313,15 +313,24 @@ size_t ASMbase::getNoElms (bool includeZeroVolElms, bool includeXElms) const
void ASMbase::getNoIntPoints (size_t& nPt, size_t& nIPt)
{
size_t nGp = 1;
for (unsigned char d = 0; d < ndim; d++)
nGp *= nGauss;
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 };
this->getOrder(ng[0],ng[1],ng[2]);
for (unsigned char d = 0; d < ndim; d++)
nGp *= ng[d] + nGauss%10;
}
firstIp = nPt;
nPt += nel*nGp; // Note: Includes also the 0-span elements
// Count additional interface quadrature points
size_t nInterface = MLGE.size() - nel;
if (nInterface > 0 && nInterface != nel && nGauss > 0)
if (nInterface > 0 && nInterface != nel && nGauss > 0 && nGauss <= 10)
nIPt += nInterface*nGp/nGauss;
}
@ -329,8 +338,18 @@ void ASMbase::getNoIntPoints (size_t& nPt, size_t& nIPt)
void ASMbase::getNoBouPoints (size_t& nPt, char ldim, char lindx)
{
size_t nGp = 1;
for (char d = 0; d < ldim; d++)
nGp *= nGauss;
if (nGauss > 0 && nGauss <= 10)
for (char d = 0; d < ldim; d++)
nGp *= nGauss;
else
{
// Use polynomial order to define number of quadrature points
int ng[3] = { 0, 0, 0 };
this->getOrder(ng[0],ng[1],ng[2]);
ng[(lindx-1)/2] = 1;
for (unsigned char d = 0; d < ndim; d++)
nGp *= ng[d];
}
firstBp[lindx] = nPt;
@ -1186,6 +1205,17 @@ int ASMbase::searchCtrlPt (RealArray::const_iterator cit,
}
int ASMbase::getNoGaussPt (int p, bool neumann) const
{
if (nGauss > 0 && nGauss < 11)
return nGauss;
else if (neumann)
return p;
return p + nGauss%10;
}
bool ASMbase::evalSolution (Matrix&, const Vector&, const int*, int) const
{
return Aerror("evalSolution(Matrix&,const Vector&,const int*,int)");

19
src/ASM/ASMbase.h
View File

@ -150,8 +150,7 @@ public:
bool addGlobalLagrangeMultipliers(const IntVec& mGLag,
unsigned char nnLag = 1);
//! \brief Defines the numerical integration scheme to use.
//! \param[in] ng Number of Gauss points in each parameter direction
//! \brief Defines the numerical integration scheme \a nGauss in the patch.
void setGauss(int ng) { nGauss = ng; }
//! \brief Defines the number of solution fields \a nf in the patch.
@ -675,6 +674,11 @@ protected:
//! \param[in] pch Pointer to the neighboring patch
void addNeighbor(ASMbase* pch);
//! \brief Returns the number of Gauss points to use in one direction.
//! \param[in] p Polynomial order of the basis functions
//! \param[in] neumann Whether or not we are assembling Neumann BCs
int getNoGaussPt(int p, bool neumann = false) const;
//! \brief Helper method used by evalPoint to search for a control point.
//! \param[in] cit iterator of array of control point coordinates
//! \param[in] end iterator of array of control point coordinates
@ -747,7 +751,6 @@ protected:
unsigned char nsd; //!< Number of space dimensions (ndim <= nsd <= 3)
unsigned char nf; //!< Number of primary solution fields (1 or larger)
unsigned char nLag; //!< Number of Lagrange multipliers per node
int nGauss; //!< Numerical integration scheme
size_t nel; //!< Number of regular elements in this patch
size_t nnod; //!< Number of regular nodes in this patch
@ -768,6 +771,16 @@ protected:
IntVec myMLGN; //!< The actual Matrix of Local to Global Node numbers
IntMat myMNPC; //!< The actual Matrix of Nodal Point Correspondance
//! \brief Numerical integration scheme for this patch.
//! \details A value in the range [1,10] means use that number of Gauss
//! quadrature points in each parameter direction, regardless of polynomial
//! order of the basis functions. If zero or negative, the number of Gauss
//! quadrature points is set independently in each parameter direction to
//! \a p+nGauss, where \a p is the polynomial order in that direction.
//! If the value is set larger than 10, the number of quadrature points
//! in each parameter direction is set to \a p+nGauss%10.
int nGauss; //!< \sa getNoGaussPt
size_t firstIp; //!< Global index to first interior integration point
//! Global indices to first integration point for the Neumann boundaries

29
src/ASM/ASMs1D.C
View File

@ -746,15 +746,18 @@ bool ASMs1D::integrate (Integrand& integrand,
{
if (!curv) return true; // silently ignore empty patches
const int p1 = curv->order();
// Get Gaussian quadrature points and weights
const double* xg = GaussQuadrature::getCoord(nGauss);
const double* wg = GaussQuadrature::getWeight(nGauss);
const int ng = this->getNoGaussPt(p1);
const double* xg = GaussQuadrature::getCoord(ng);
const double* wg = GaussQuadrature::getWeight(ng);
if (!xg || !wg) return false;
// Get the reduced integration quadrature points, if needed
const double* xr = nullptr;
const double* wr = nullptr;
int nRed = integrand.getReducedIntegration(nGauss);
int nRed = integrand.getReducedIntegration(ng);
if (nRed > 0)
{
xr = GaussQuadrature::getCoord(nRed);
@ -762,7 +765,7 @@ bool ASMs1D::integrate (Integrand& integrand,
if (!xr || !wr) return false;
}
else if (nRed < 0)
nRed = nGauss; // The integrand needs to know nGauss
nRed = ng; // The integrand needs to know nGauss
if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES)
if (curv->rational())
@ -774,12 +777,10 @@ bool ASMs1D::integrate (Integrand& integrand,
// Compute parameter values of the Gauss points over the whole patch
Matrix gpar, redpar;
this->getGaussPointParameters(gpar,nGauss,xg);
this->getGaussPointParameters(gpar,ng,xg);
if (xr)
this->getGaussPointParameters(redpar,nRed,xr);
const int p1 = curv->order();
FiniteElement fe(p1);
Matrix dNdu, Jac;
Matrix3D d2Ndu2, Hess;
@ -852,10 +853,10 @@ bool ASMs1D::integrate (Integrand& integrand,
// --- Integration loop over all Gauss points in current element -----------
int jp = iel*nGauss;
int jp = iel*ng;
fe.iGP = firstIp + jp; // Global integration point counter
for (int i = 0; i < nGauss && ok; i++, fe.iGP++)
for (int i = 0; i < ng && ok; i++, fe.iGP++)
{
// Local element coordinate of current integration point
fe.xi = xg[i];
@ -1368,16 +1369,16 @@ bool ASMs1D::assembleL2matrices (SparseMatrix& A, StdVector& B,
const int p1 = curv->order();
// Get Gaussian quadrature point coordinates (and weights if continuous)
const int ng1 = continuous ? nGauss : p1 - 1;
const double* xg = GaussQuadrature::getCoord(ng1);
const double* wg = continuous ? GaussQuadrature::getWeight(nGauss) : nullptr;
const int ng = continuous ? this->getNoGaussPt(p1,true) : p1 - 1;
const double* xg = GaussQuadrature::getCoord(ng);
const double* wg = continuous ? GaussQuadrature::getWeight(ng) : nullptr;
if (!xg) return false;
if (continuous && !wg) return false;
// Compute parameter values of the Gauss points over the whole patch
// and evaluate the secondary solution at all integration points
Matrix gp, sField;
RealArray gpar = this->getGaussPointParameters(gp,ng1,xg);
RealArray gpar = this->getGaussPointParameters(gp,ng,xg);
if (!this->evalSolution(sField,integrand,&gpar))
{
std::cerr <<" *** ASMs1D::assembleL2matrices: Failed for patch "<< idx+1
@ -1408,7 +1409,7 @@ bool ASMs1D::assembleL2matrices (SparseMatrix& A, StdVector& B,
// --- Integration loop over all Gauss points in current element -----------
for (int i = 0; i < ng1; i++, ip++)
for (int i = 0; i < ng; i++, ip++)
{
// Fetch basis function values at current integration point
if (continuous)

13
src/ASM/ASMs1DLag.C
View File

@ -202,14 +202,15 @@ bool ASMs1DLag::integrate (Integrand& integrand,
if (this->empty()) return true; // silently ignore empty patches
// Get Gaussian quadrature points and weights
const double* xg = GaussQuadrature::getCoord(nGauss);
const double* wg = GaussQuadrature::getWeight(nGauss);
const int ng = this->getNoGaussPt(p1);
const double* xg = GaussQuadrature::getCoord(ng);
const double* wg = GaussQuadrature::getWeight(ng);
if (!xg || !wg) return false;
// Get the reduced integration quadrature points, if needed
const double* xr = nullptr;
const double* wr = nullptr;
int nRed = integrand.getReducedIntegration(nGauss);
int nRed = integrand.getReducedIntegration(ng);
if (nRed > 0)
{
xr = GaussQuadrature::getCoord(nRed);
@ -217,7 +218,7 @@ bool ASMs1DLag::integrate (Integrand& integrand,
if (!xr || !wr) return false;
}
else if (nRed < 0)
nRed = nGauss; // The integrand needs to know nGauss
nRed = ng; // The integrand needs to know nGauss
// Get parametric coordinates of the elements
RealArray gpar;
@ -285,10 +286,10 @@ bool ASMs1DLag::integrate (Integrand& integrand,
// --- Integration loop over all Gauss points in each direction ------------
int jp = iel*nGauss;
int jp = iel*ng;
fe.iGP = firstIp + jp; // Global integration point counter
for (int i = 0; i < nGauss && ok; i++, fe.iGP++)
for (int i = 0; i < ng && ok; i++, fe.iGP++)
{
// Local element coordinate of current integration point
fe.xi = xg[i];

93
src/ASM/ASMs2D.C
View File

@ -1486,15 +1486,24 @@ bool ASMs2D::integrate (Integrand& integrand,
bool use2ndDer = integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES;
bool useElmVtx = integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS;
const int p1 = surf->order_u();
const int p2 = surf->order_v();
// Get Gaussian quadrature points and weights
const double* xg = GaussQuadrature::getCoord(nGauss);
const double* wg = GaussQuadrature::getWeight(nGauss);
if (!xg || !wg) return false;
std::array<int,2> ng;
std::array<const double*,2> xg, wg;
for (int d = 0; d < 2; d++)
{
ng[d] = this->getNoGaussPt(d == 0 ? p1 : p2);
xg[d] = GaussQuadrature::getCoord(ng[d]);
wg[d] = GaussQuadrature::getWeight(ng[d]);
if (!xg[d] || !wg[d]) return false;
}
// Get the reduced integration quadrature points, if needed
const double* xr = nullptr;
const double* wr = nullptr;
int nRed = integrand.getReducedIntegration(nGauss);
int nRed = integrand.getReducedIntegration(ng[0]);
if (nRed > 0)
{
xr = GaussQuadrature::getCoord(nRed);
@ -1502,13 +1511,13 @@ bool ASMs2D::integrate (Integrand& integrand,
if (!xr && !wr) return false;
}
else if (nRed < 0)
nRed = nGauss; // The integrand needs to know nGauss
nRed = ng[0]; // The integrand needs to know nGauss
// Compute parameter values of the Gauss points over the whole patch
std::array<Matrix,2> gpar, redpar;
for (int d = 0; d < 2; d++)
{
this->getGaussPointParameters(gpar[d],d,nGauss,xg);
this->getGaussPointParameters(gpar[d],d,ng[d],xg[d]);
if (xr)
this->getGaussPointParameters(redpar[d],d,nRed,xr);
}
@ -1529,8 +1538,6 @@ bool ASMs2D::integrate (Integrand& integrand,
std::cout <<"\nBasis functions at integration point "<< 1+i << spline[i];
#endif
const int p1 = surf->order_u();
const int p2 = surf->order_v();
const int n1 = surf->numCoefs_u();
const int nel1 = n1 - p1 + 1;
@ -1596,9 +1603,9 @@ bool ASMs2D::integrate (Integrand& integrand,
fe.Navg.resize(p1*p2,true);
double area = 0.0;
int ip = ((i2-p2)*nGauss*nel1 + i1-p1)*nGauss;
for (int j = 0; j < nGauss; j++, ip += nGauss*(nel1-1))
for (int i = 0; i < nGauss; i++, ip++)
int ip = ((i2-p2)*ng[1]*nel1 + i1-p1)*ng[0];
for (int j = 0; j < ng[1]; j++, ip += ng[0]*(nel1-1))
for (int i = 0; i < ng[0]; i++, ip++)
{
// Fetch basis function derivatives at current integration point
SplineUtils::extractBasis(spline[ip],fe.N,dNdu);
@ -1606,7 +1613,7 @@ bool ASMs2D::integrate (Integrand& integrand,
// Compute Jacobian determinant of coordinate mapping
// and multiply by weight of current integration point
double detJac = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu,false);
double weight = dA*wg[i]*wg[j];
double weight = dA*wg[0][i]*wg[1][j];
// Numerical quadrature
fe.Navg.add(fe.N,detJac*weight);
@ -1620,8 +1627,8 @@ bool ASMs2D::integrate (Integrand& integrand,
if (integrand.getIntegrandType() & Integrand::ELEMENT_CENTER)
{
// Compute the element center
double u0 = 0.5*(gpar[0](1,i1-p1+1) + gpar[0](nGauss,i1-p1+1));
double v0 = 0.5*(gpar[1](1,i2-p2+1) + gpar[1](nGauss,i2-p2+1));
double u0 = 0.5*(gpar[0](1,i1-p1+1) + gpar[0](ng[0],i1-p1+1));
double v0 = 0.5*(gpar[1](1,i2-p2+1) + gpar[1](ng[1],i2-p2+1));
SplineUtils::point(X,u0,v0,surf);
if (!useElmVtx)
{
@ -1682,16 +1689,16 @@ bool ASMs2D::integrate (Integrand& integrand,
// --- Integration loop over all Gauss points in each direction --------
int ip = ((i2-p2)*nGauss*nel1 + i1-p1)*nGauss;
int jp = ((i2-p2)*nel1 + i1-p1)*nGauss*nGauss;
int ip = ((i2-p2)*ng[1]*nel1 + i1-p1)*ng[0];
int jp = ((i2-p2)*nel1 + i1-p1)*ng[0]*ng[1];
fe.iGP = firstIp + jp; // Global integration point counter
for (int j = 0; j < nGauss; j++, ip += nGauss*(nel1-1))
for (int i = 0; i < nGauss; i++, ip++, fe.iGP++)
for (int j = 0; j < ng[1]; j++, ip += ng[0]*(nel1-1))
for (int i = 0; i < ng[0]; i++, ip++, fe.iGP++)
{
// Local element coordinates of current integration point
fe.xi = xg[i];
fe.eta = xg[j];
fe.xi = xg[0][i];
fe.eta = xg[1][j];
// Parameter values of current integration point
fe.u = param[0] = gpar[0](i+1,i1-p1+1);
@ -1724,16 +1731,7 @@ bool ASMs2D::integrate (Integrand& integrand,
#if SP_DEBUG > 4
if (iel == dbgElm || iel == -dbgElm || dbgElm == 0)
{
std::cout <<"\niel, ip = "<< iel <<" "<< ip
<<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
if (!fe.d2NdX2.empty())
std::cout <<"d2NdX2 ="<< fe.d2NdX2;
if (!fe.G.empty())
std::cout <<"G ="<< fe.G;
if (!fe.H.empty())
std::cout <<"H ="<< fe.H;
}
std::cout <<"\n"<< fe;
#endif
// Cartesian coordinates of current integration point
@ -1741,7 +1739,7 @@ bool ASMs2D::integrate (Integrand& integrand,
X.t = time.t;
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= dA*wg[i]*wg[j];
fe.detJxW *= dA*wg[0][i]*wg[1][j];
#ifndef USE_OPENMP
PROFILE3("Integrand::evalInt");
#endif
@ -1777,7 +1775,7 @@ bool ASMs2D::integrate (Integrand& integrand,
{
if (!surf) return true; // silently ignore empty patches
if (integrand.getReducedIntegration(nGauss) != 0)
if (integrand.getReducedIntegration(2) != 0)
{
std::cerr <<" *** ASMs2D::integrate(Integrand&,GlobalIntegral&,"
<<"const TimeDomain&,const Real3DMat&): Available for standard"
@ -1926,8 +1924,7 @@ bool ASMs2D::integrate (Integrand& integrand,
#if SP_DEBUG > 4
if (iel == dbgElm || iel == -dbgElm || dbgElm == 0)
std::cout <<"\niel, jp = "<< iel <<" "<< jp
<<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
std::cout <<"\n"<< fe;
#endif
// Cartesian coordinates of current integration point
@ -1975,13 +1972,15 @@ bool ASMs2D::integrate (Integrand& integrand,
PROFILE2("ASMs2D::integrate(J)");
// Get Gaussian quadrature points and weights
const double* xg = GaussQuadrature::getCoord(nGauss);
const double* wg = GaussQuadrature::getWeight(nGauss);
if (!xg || !wg) return false;
const int p1 = surf->order_u();
const int p2 = surf->order_v();
const int ng = this->getNoGaussPt(p1 > p2 ? p1 : p2);
// Get Gaussian quadrature points and weights
const double* xg = GaussQuadrature::getCoord(ng);
const double* wg = GaussQuadrature::getWeight(ng);
if (!xg || !wg) return false;
const int n1 = surf->numCoefs_u();
const int n2 = surf->numCoefs_v();
@ -2049,7 +2048,7 @@ bool ASMs2D::integrate (Integrand& integrand,
// --- Integration loop over all Gauss points along the edge ---------
for (int i = 0; i < nGauss && ok; i++)
for (int i = 0; i < ng && ok; i++)
{
// Local element coordinates and parameter values
// of current integration point
@ -2100,13 +2099,11 @@ bool ASMs2D::integrate (Integrand& integrand,
this->extractBasis(fe.u,fe.v,t1,fe.p, dN, edgeDir > 0);
utl::merge(fe.N,dN,MNPC[iel],MNPC[kel]);
}
}
#if SP_DEBUG > 4
std::cout <<"\niel, xi,eta = "<< fe.iel
<<" "<< fe.xi <<" "<< fe.eta
<<"\ndN ="<< fe.N <<"dNdX ="<< fe.dNdX;
std::cout <<"\n"<< fe;
#endif
}
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= dS*wg[i];
@ -2135,8 +2132,12 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex,
PROFILE2("ASMs2D::integrate(B)");
const int p1 = surf->order_u();
const int p2 = surf->order_v();
// Get Gaussian quadrature points and weights
int nGP = integrand.getBouIntegrationPoints(nGauss);
int nG1 = this->getNoGaussPt(lIndex%10 < 3 ? p1 : p2, true);
int nGP = integrand.getBouIntegrationPoints(nG1);
const double* xg = GaussQuadrature::getCoord(nGP);
const double* wg = GaussQuadrature::getWeight(nGP);
if (!xg || !wg) return false;
@ -2170,8 +2171,6 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex,
std::vector<Go::BasisDerivsSf> spline;
surf->computeBasisGrid(gpar[0],gpar[1],spline);
const int p1 = surf->order_u();
const int p2 = surf->order_v();
const int n1 = surf->numCoefs_u();
const int n2 = surf->numCoefs_v();

38
src/ASM/ASMs2DLag.C
View File

@ -25,6 +25,7 @@
#include "ElementBlock.h"
#include "Utilities.h"
#include "Vec3Oper.h"
#include <array>
ASMs2DLag::ASMs2DLag (unsigned char n_s, unsigned char n_f)
@ -306,14 +307,20 @@ bool ASMs2DLag::integrate (Integrand& integrand,
if (this->empty()) return true; // silently ignore empty patches
// Get Gaussian quadrature points and weights
const double* xg = GaussQuadrature::getCoord(nGauss);
const double* wg = GaussQuadrature::getWeight(nGauss);
if (!xg || !wg) return false;
std::array<int,2> ng;
std::array<const double*,2> xg, wg;
for (int d = 0; d < 2; d++)
{
ng[d] = this->getNoGaussPt(d == 0 ? p1 : p2);
xg[d] = GaussQuadrature::getCoord(ng[d]);
wg[d] = GaussQuadrature::getWeight(ng[d]);
if (!xg[d] || !wg[d]) return false;
}
// Get the reduced integration quadrature points, if needed
const double* xr = nullptr;
const double* wr = nullptr;
int nRed = integrand.getReducedIntegration(nGauss);
int nRed = integrand.getReducedIntegration(ng[0]);
if (nRed > 0)
{
xr = GaussQuadrature::getCoord(nRed);
@ -321,7 +328,7 @@ bool ASMs2DLag::integrate (Integrand& integrand,
if (!xr || !wr) return false;
}
else if (nRed < 0)
nRed = nGauss; // The integrand needs to know nGauss
nRed = ng[0]; // The integrand needs to know nGauss
// Get parametric coordinates of the elements
RealArray upar, vpar;
@ -416,25 +423,25 @@ bool ASMs2DLag::integrate (Integrand& integrand,
// --- Integration loop over all Gauss points in each direction --------
int jp = iel*nGauss*nGauss;
int jp = iel*ng[0]*ng[1];
fe.iGP = firstIp + jp; // Global integration point counter
for (int j = 0; j < nGauss; j++)
for (int i = 0; i < nGauss; i++, fe.iGP++)
for (int j = 0; j < ng[1]; j++)
for (int i = 0; i < ng[0]; i++, fe.iGP++)
{
// Local element coordinates of current integration point
fe.xi = xg[i];
fe.eta = xg[j];
fe.xi = xg[0][i];
fe.eta = xg[1][j];
// Parameter value of current integration point
if (!upar.empty())
fe.u = 0.5*(upar[i1]*(1.0-xg[i]) + upar[i1+1]*(1.0+xg[i]));
fe.u = 0.5*(upar[i1]*(1.0-xg[0][i]) + upar[i1+1]*(1.0+xg[0][i]));
if (!vpar.empty())
fe.v = 0.5*(vpar[i2]*(1.0-xg[j]) + vpar[i2+1]*(1.0+xg[j]));
fe.v = 0.5*(vpar[i2]*(1.0-xg[1][j]) + vpar[i2+1]*(1.0+xg[1][j]));
// Compute basis function derivatives at current integration point
// using tensor product of one-dimensional Lagrange polynomials
if (!Lagrange::computeBasis(fe.N,dNdu,p1,xg[i],p2,xg[j]))
if (!Lagrange::computeBasis(fe.N,dNdu,p1,xg[0][i],p2,xg[1][j]))
ok = false;
// Compute Jacobian inverse of coordinate mapping and derivatives
@ -446,7 +453,7 @@ bool ASMs2DLag::integrate (Integrand& integrand,
X.t = time.t;
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg[i]*wg[j];
fe.detJxW *= wg[0][i]*wg[1][j];
if (!integrand.evalInt(*A,fe,time,X))
ok = false;
}
@ -475,7 +482,8 @@ bool ASMs2DLag::integrate (Integrand& integrand, int lIndex,
if (this->empty()) return true; // silently ignore empty patches
// Get Gaussian quadrature points and weights
int nGP = integrand.getBouIntegrationPoints(nGauss);
int nG1 = this->getNoGaussPt(lIndex%10 < 3 ? p1 : p2, true);
int nGP = integrand.getBouIntegrationPoints(nG1);
const double* xg = GaussQuadrature::getCoord(nGP);
const double* wg = GaussQuadrature::getWeight(nGP);
if (!xg || !wg) return false;

21
src/ASM/ASMs2Drecovery.C
View File

@ -172,13 +172,14 @@ bool ASMs2D::assembleL2matrices (SparseMatrix& A, StdVector& B,
const int n2 = surf->numCoefs_v();
const int nel1 = n1 - p1 + 1;
const int nel2 = n2 - p2 + 1;
const int pmax = p1 > p2 ? p1 : p2;
// Get Gaussian quadrature point coordinates (and weights if continuous)
const int ng1 = continuous ? nGauss : p1 - 1;
const int ng2 = continuous ? nGauss : p2 - 1;
const int ng1 = continuous ? this->getNoGaussPt(pmax,true) : p1 - 1;
const int ng2 = continuous ? ng1 : p2 - 1;
const double* xg = GaussQuadrature::getCoord(ng1);
const double* yg = GaussQuadrature::getCoord(ng2);
const double* wg = continuous ? GaussQuadrature::getWeight(nGauss) : nullptr;
const double* wg = continuous ? GaussQuadrature::getWeight(ng1) : nullptr;
if (!xg || !yg) return false;
if (continuous && !wg) return false;
@ -505,8 +506,12 @@ Go::SplineSurface* ASMs2D::projectSolutionLeastSquare (const IntegrandBase& inte
if (!surf) return nullptr;
// Get Gaussian quadrature points and weights
const double* xg = GaussQuadrature::getCoord(nGauss);
const double* wg = GaussQuadrature::getWeight(nGauss);
const int p1 = surf->order_u();
const int p2 = surf->order_v();
const int ng = this->getNoGaussPt(p1 > p2 ? p1 : p2, true);
const double* xg = GaussQuadrature::getCoord(ng);
const double* wg = GaussQuadrature::getWeight(ng);
if (!xg || !wg) return nullptr;
// Compute parameter values of the result sampling points (Gauss points)
@ -514,18 +519,18 @@ Go::SplineSurface* ASMs2D::projectSolutionLeastSquare (const IntegrandBase& inte
std::array<RealArray,2> gpar, wgpar;
for (int dir = 0; dir < 2; dir++)
{
this->getGaussPointParameters(ggpar[dir],dir,nGauss,xg);
this->getGaussPointParameters(ggpar[dir],dir,ng,xg);
gpar[dir] = ggpar[dir];
// Gauss weights at parameter values
const Go::BsplineBasis& basis = surf->basis(dir);
RealArray::const_iterator knotit = basis.begin();
RealArray& tmp = wgpar[dir];
tmp.reserve(nGauss*(basis.numCoefs()-basis.order()));
tmp.reserve(ng*(basis.numCoefs()-basis.order()));
for (int i = basis.order(); i <= basis.numCoefs(); i++)
{
double d = knotit[i] - knotit[i-1];
for (int j = 0; j < nGauss; j++)
for (int j = 0; j < ng; j++)
tmp.push_back(d > 0.0 ? wg[j]*d*0.5 : 0.0);
}
}

97
src/ASM/ASMs3D.C
View File

@ -1677,14 +1677,20 @@ bool ASMs3D::integrate (Integrand& integrand,
bool useElmVtx = integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS;
// Get Gaussian quadrature points and weights
const double* xg = GaussQuadrature::getCoord(nGauss);
const double* wg = GaussQuadrature::getWeight(nGauss);
if (!xg || !wg) return false;
std::array<int,3> ng;
std::array<const double*,3> xg, wg;
for (int d = 0; d < 3; d++)
{
ng[d] = this->getNoGaussPt(svol->order(d));
xg[d] = GaussQuadrature::getCoord(ng[d]);
wg[d] = GaussQuadrature::getWeight(ng[d]);
if (!xg[d] || !wg[d]) return false;
}
// Get the reduced integration quadrature points, if needed
const double* xr = nullptr;
const double* wr = nullptr;
int nRed = integrand.getReducedIntegration(nGauss);
int nRed = integrand.getReducedIntegration(ng[0]);
if (nRed > 0)
{
xr = GaussQuadrature::getCoord(nRed);
@ -1692,13 +1698,13 @@ bool ASMs3D::integrate (Integrand& integrand,
if (!xr || !wr) return false;
}
else if (nRed < 0)
nRed = nGauss; // The integrand needs to know nGauss
nRed = ng[0]; // The integrand needs to know nGauss
// Compute parameter values of the Gauss points over the whole patch
std::array<Matrix,3> gpar, redpar;
for (int d = 0; d < 3; d++)
{
this->getGaussPointParameters(gpar[d],d,nGauss,xg);
this->getGaussPointParameters(gpar[d],d,ng[d],xg[d]);
if (xr)
this->getGaussPointParameters(redpar[d],d,nRed,xr);
}
@ -1789,10 +1795,10 @@ bool ASMs3D::integrate (Integrand& integrand,
fe.Navg.resize(p1*p2*p3,true);
double vol = 0.0;
int ip = (((i3-p3)*nGauss*nel2 + i2-p2)*nGauss*nel1 + i1-p1)*nGauss;
for (int k = 0; k < nGauss; k++, ip += nGauss*(nel2-1)*nGauss*nel1)
for (int j = 0; j < nGauss; j++, ip += nGauss*(nel1-1))
for (int i = 0; i < nGauss; i++, ip++)
int ip = (((i3-p3)*ng[2]*nel2 + i2-p2)*ng[1]*nel1 + i1-p1)*ng[0];
for (int k = 0; k < ng[2]; k++, ip += ng[1]*(nel2-1)*ng[0]*nel1)
for (int j = 0; j < ng[1]; j++, ip += ng[0]*(nel1-1))
for (int i = 0; i < ng[0]; i++, ip++)
{
// Fetch basis function derivatives at current integration point
SplineUtils::extractBasis(spline[ip],fe.N,dNdu);
@ -1800,7 +1806,7 @@ bool ASMs3D::integrate (Integrand& integrand,
// Compute Jacobian determinant of coordinate mapping
// and multiply by weight of current integration point
double detJac = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu,false);
double weight = 0.125*dV*wg[i]*wg[j]*wg[k];
double weight = 0.125*dV*wg[0][i]*wg[1][j]*wg[2][k];
// Numerical quadrature
fe.Navg.add(fe.N,detJac*weight);
@ -1814,9 +1820,9 @@ bool ASMs3D::integrate (Integrand& integrand,
if (integrand.getIntegrandType() & Integrand::ELEMENT_CENTER)
{
// Compute the element center
double u0 = 0.5*(gpar[0](1,i1-p1+1) + gpar[0](nGauss,i1-p1+1));
double v0 = 0.5*(gpar[1](1,i2-p2+1) + gpar[1](nGauss,i2-p2+1));
double w0 = 0.5*(gpar[2](1,i3-p3+1) + gpar[2](nGauss,i3-p3+1));
double u0 = 0.5*(gpar[0](1,i1-p1+1) + gpar[0](ng[0],i1-p1+1));
double v0 = 0.5*(gpar[1](1,i2-p2+1) + gpar[1](ng[1],i2-p2+1));
double w0 = 0.5*(gpar[2](1,i3-p3+1) + gpar[2](ng[2],i3-p3+1));
SplineUtils::point(X,u0,v0,w0,svol);
if (!useElmVtx)
{
@ -1876,18 +1882,18 @@ bool ASMs3D::integrate (Integrand& integrand,
// --- Integration loop over all Gauss points in each direction --------
int ip = (((i3-p3)*nGauss*nel2 + i2-p2)*nGauss*nel1 + i1-p1)*nGauss;
int jp = (((i3-p3)*nel2 + i2-p2)*nel1 + i1-p1)*nGauss*nGauss*nGauss;
int ip = (((i3-p3)*ng[2]*nel2 + i2-p2)*ng[1]*nel1 + i1-p1)*ng[0];
int jp = (((i3-p3)*nel2 + i2-p2)*nel1 + i1-p1)*ng[0]*ng[1]*ng[2];
fe.iGP = firstIp + jp; // Global integration point counter
for (int k = 0; k < nGauss; k++, ip += nGauss*(nel2-1)*nGauss*nel1)
for (int j = 0; j < nGauss; j++, ip += nGauss*(nel1-1))
for (int i = 0; i < nGauss; i++, ip++, fe.iGP++)
for (int k = 0; k < ng[2]; k++, ip += ng[1]*(nel2-1)*ng[0]*nel1)
for (int j = 0; j < ng[1]; j++, ip += ng[0]*(nel1-1))
for (int i = 0; i < ng[0]; i++, ip++, fe.iGP++)
{
// Local element coordinates of current integration point
fe.xi = xg[i];
fe.eta = xg[j];
fe.zeta = xg[k];
fe.xi = xg[0][i];
fe.eta = xg[1][j];
fe.zeta = xg[2][k];
// Parameter values of current integration point
fe.u = param[0] = gpar[0](i+1,i1-p1+1);
@ -1915,12 +1921,7 @@ bool ASMs3D::integrate (Integrand& integrand,
#if SP_DEBUG > 4
if (iel == dbgElm || iel == -dbgElm || dbgElm == 0)
{
std::cout <<"\niel, ip = "<< iel <<" "<< ip
<<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
if (!fe.d2NdX2.empty())
std::cout <<"d2NdX2 ="<< fe.d2NdX2;
}
std::cout <<"\n"<< fe;
#endif
// Cartesian coordinates of current integration point
@ -1928,7 +1929,7 @@ bool ASMs3D::integrate (Integrand& integrand,
X.t = time.t;
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= 0.125*dV*wg[i]*wg[j]*wg[k];
fe.detJxW *= 0.125*dV*wg[0][i]*wg[1][j]*wg[2][k];
#ifndef USE_OPENMP
PROFILE3("Integrand::evalInt");
#endif
@ -1964,7 +1965,7 @@ bool ASMs3D::integrate (Integrand& integrand,
{
if (!svol) return true; // silently ignore empty patches
if (integrand.getReducedIntegration(nGauss) != 0)
if (integrand.getReducedIntegration(2) != 0)
{
std::cerr <<" *** ASMs3D::integrate(Integrand&,GlobalIntegral&,"
<<"const TimeDomain&,const Real3DMat&): Available for standard"
@ -2112,8 +2113,7 @@ bool ASMs3D::integrate (Integrand& integrand,
#if SP_DEBUG > 4
if (iel == dbgElm || iel == -dbgElm || dbgElm == 0)
std::cout <<"\niel, jp = "<< iel <<" "<< jp
<<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
std::cout <<"\n"<< fe;
#endif
// Cartesian coordinates of current integration point
@ -2181,18 +2181,20 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex,
}
const ThreadGroups& threadGrp = tit->second;
// Get Gaussian quadrature points and weights
int nGP = integrand.getBouIntegrationPoints(nGauss);
const double* xg = GaussQuadrature::getCoord(nGP);
const double* wg = GaussQuadrature::getWeight(nGP);
if (!xg || !wg) return false;
// Find the parametric direction of the face normal {-3,-2,-1, 1, 2, 3}
const int faceDir = (lIndex%10+1)/(lIndex%2 ? -2 : 2);
const int t1 = 1 + abs(faceDir)%3; // first tangent direction
const int t2 = 1 + t1%3; // second tangent direction
// Get Gaussian quadrature points and weights
// For now, use the largest polynomial order of the two tangent directions
int nG1 = this->getNoGaussPt(std::max(svol->order(t1),svol->order(t2)),true);
int nGP = integrand.getBouIntegrationPoints(nG1);
const double* xg = GaussQuadrature::getCoord(nGP);
const double* wg = GaussQuadrature::getWeight(nGP);
if (!xg || !wg) return false;
// Compute parameter values of the Gauss points over the whole patch face
std::array<Matrix,3> gpar;
for (int d = 0; d < 3; d++)
@ -2414,8 +2416,9 @@ bool ASMs3D::integrateEdge (Integrand& integrand, int lEdge,
PROFILE2("ASMs3D::integrate(E)");
// Get Gaussian quadrature points and weights
const double* xg = GaussQuadrature::getCoord(nGauss);
const double* wg = GaussQuadrature::getWeight(nGauss);
int ng = this->getNoGaussPt(svol->order((lEdge-1)/4),true);
const double* xg = GaussQuadrature::getCoord(ng);
const double* wg = GaussQuadrature::getWeight(ng);
if (!xg || !wg) return false;
// Compute parameter values of the Gauss points along the whole patch edge
@ -2441,11 +2444,11 @@ bool ASMs3D::integrateEdge (Integrand& integrand, int lEdge,
RealArray::const_iterator uit = svol->basis(d).begin() + pm1;
double ucurr, uprev = *(uit++);
int nCol = svol->numCoefs(d) - pm1;
gpar[d].resize(nGauss,nCol);
gpar[d].resize(ng,nCol);
for (int j = 1; j <= nCol; ++uit, j++)
{
ucurr = *uit;
for (int i = 1; i <= nGauss; i++)
for (int i = 1; i <= ng; i++)
gpar[d](i,j) = 0.5*((ucurr-uprev)*xg[i-1] + ucurr+uprev);
uprev = ucurr;
}
@ -2521,17 +2524,17 @@ bool ASMs3D::integrateEdge (Integrand& integrand, int lEdge,
if (lEdge < 5)
{
dS = svol->knotSpan(0,nodeInd[ip].I);
ip = (i1-p1)*nGauss;
ip = (i1-p1)*ng;
}
else if (lEdge < 9)
{
dS = svol->knotSpan(1,nodeInd[ip].J);
ip = (i2-p2)*nGauss;
ip = (i2-p2)*ng;
}
else if (lEdge < 13)
{
dS = svol->knotSpan(2,nodeInd[ip].K);
ip = (i3-p3)*nGauss;
ip = (i3-p3)*ng;
}
// Set up control point coordinates for current element
@ -2546,7 +2549,7 @@ bool ASMs3D::integrateEdge (Integrand& integrand, int lEdge,
fe.iGP = firstp + ip; // Global integration point counter
for (int i = 0; i < nGauss && ok; i++, ip++, fe.iGP++)
for (int i = 0; i < ng && ok; i++, ip++, fe.iGP++)
{
// Parameter values of current integration point
if (gpar[0].size() > 1) fe.u = param[0] = gpar[0](i+1,i1-p1+1);
@ -2556,7 +2559,7 @@ bool ASMs3D::integrateEdge (Integrand& integrand, int lEdge,
// Fetch basis function derivatives at current integration point
SplineUtils::extractBasis(spline[ip],fe.N,dNdu);
// Compute basis function derivatives and the edge tang
// Compute basis function derivatives and the edge tangent
fe.detJxW = utl::Jacobian(Jac,tang,fe.dNdX,Xnod,dNdu,1+(lEdge-1)/4);
if (fe.detJxW == 0.0) continue; // skip singular points

71
src/ASM/ASMs3DLag.C
View File

@ -337,14 +337,20 @@ bool ASMs3DLag::integrate (Integrand& integrand,
if (this->empty()) return true; // silently ignore empty patches
// Get Gaussian quadrature points and weights
const double* xg = GaussQuadrature::getCoord(nGauss);
const double* wg = GaussQuadrature::getWeight(nGauss);
if (!xg || !wg) return false;
std::array<int,3> ng;
std::array<const double*,3> xg, wg;
for (int d = 0; d < 3; d++)
{
ng[d] = this->getNoGaussPt(d == 0 ? p1 : (d == 1 ? p2 : p3));
xg[d] = GaussQuadrature::getCoord(ng[d]);
wg[d] = GaussQuadrature::getWeight(ng[d]);
if (!xg[d] || !wg[d]) return false;
}
// Get the reduced integration quadrature points, if needed
const double* xr = nullptr;
const double* wr = nullptr;
int nRed = integrand.getReducedIntegration(nGauss);
int nRed = integrand.getReducedIntegration(ng[0]);
if (nRed > 0)
{
xr = GaussQuadrature::getCoord(nRed);
@ -352,7 +358,7 @@ bool ASMs3DLag::integrate (Integrand& integrand,
if (!xr || !wr) return false;
}
else if (nRed < 0)
nRed = nGauss; // The integrand needs to know nGauss
nRed = ng[0]; // The integrand needs to know nGauss
// Get parametric coordinates of the elements
RealArray upar, vpar, wpar;
@ -458,29 +464,30 @@ bool ASMs3DLag::integrate (Integrand& integrand,
// --- Integration loop over all Gauss points in each direction --------
int jp = iel*nGauss*nGauss*nGauss;
int jp = iel*ng[0]*ng[1]*ng[2];
fe.iGP = firstIp + jp; // Global integration point counter
for (int k = 0; k < nGauss; k++)
for (int j = 0; j < nGauss; j++)
for (int i = 0; i < nGauss; i++, fe.iGP++)
for (int k = 0; k < ng[2]; k++)
for (int j = 0; j < ng[1]; j++)
for (int i = 0; i < ng[0]; i++, fe.iGP++)
{
// Local element coordinates of current integration point
fe.xi = xg[i];
fe.eta = xg[j];
fe.zeta = xg[k];
fe.xi = xg[0][i];
fe.eta = xg[1][j];
fe.zeta = xg[2][k];
// Parameter value of current integration point
if (!upar.empty())
fe.u = 0.5*(upar[i1]*(1.0-xg[i]) + upar[i1+1]*(1.0+xg[i]));
fe.u = 0.5*(upar[i1]*(1.0-fe.xi) + upar[i1+1]*(1.0+fe.xi));
if (!vpar.empty())
fe.v = 0.5*(vpar[i2]*(1.0-xg[j]) + vpar[i2+1]*(1.0+xg[j]));
fe.v = 0.5*(vpar[i2]*(1.0-fe.eta) + vpar[i2+1]*(1.0+fe.eta));
if (!wpar.empty())
fe.w = 0.5*(wpar[i3]*(1.0-xg[k]) + wpar[i3+1]*(1.0+xg[k]));
fe.w = 0.5*(wpar[i3]*(1.0-fe.zeta) + wpar[i3+1]*(1.0+fe.zeta));
// Compute basis function derivatives at current integration point
// using tensor product of one-dimensional Lagrange polynomials
if (!Lagrange::computeBasis(fe.N,dNdu,p1,xg[i],p2,xg[j],p3,xg[k]))
if (!Lagrange::computeBasis(fe.N,dNdu,
p1,fe.xi,p2,fe.eta,p3,fe.zeta))
ok = false;
// Compute Jacobian inverse of coordinate mapping and derivatives
@ -492,7 +499,7 @@ bool ASMs3DLag::integrate (Integrand& integrand,
X.t = time.t;
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg[i]*wg[j]*wg[k];
fe.detJxW *= wg[0][i]*wg[1][j]*wg[2][k];
if (!integrand.evalInt(*A,fe,time,X))
ok = false;
}
@ -529,12 +536,6 @@ bool ASMs3DLag::integrate (Integrand& integrand, int lIndex,
}
const ThreadGroups& threadGrp = tit->second;
// Get Gaussian quadrature points and weights
int nGP = integrand.getBouIntegrationPoints(nGauss);
const double* xg = GaussQuadrature::getCoord(nGP);
const double* wg = GaussQuadrature::getWeight(nGP);
if (!xg || !wg) return false;
// Find the parametric direction of the face normal {-3,-2,-1, 1, 2, 3}
const int faceDir = (lIndex%10+1)/(lIndex%2 ? -2 : 2);
@ -542,6 +543,14 @@ bool ASMs3DLag::integrate (Integrand& integrand, int lIndex,
const int t1 = 1 + t0%3; // first tangent direction of the face
const int t2 = 1 + t1%3; // second tangent direction of the face
// Get Gaussian quadrature points and weights
// For now, use the largest polynomial order of the two tangent directions
int nG1 = this->getNoGaussPt(std::max(svol->order(t1),svol->order(t2)),true);
int nGP = integrand.getBouIntegrationPoints(nG1);
const double* xg = GaussQuadrature::getCoord(nGP);
const double* wg = GaussQuadrature::getWeight(nGP);
if (!xg || !wg) return false;
// Number of elements in each direction
const int nel1 = (nx-1)/(p1-1);
const int nel2 = (ny-1)/(p2-1);
@ -711,8 +720,9 @@ bool ASMs3DLag::integrateEdge (Integrand& integrand, int lEdge,
const int lDir = (lEdge-1)/4;
// Get Gaussian quadrature points and weights
const double* xg = GaussQuadrature::getCoord(nGauss);
const double* wg = GaussQuadrature::getWeight(nGauss);
int ng = this->getNoGaussPt(svol->order((lEdge-1)/4),true);
const double* xg = GaussQuadrature::getCoord(ng);
const double* wg = GaussQuadrature::getWeight(ng);
if (!xg || !wg) return false;
// Number of elements in each direction
@ -773,11 +783,11 @@ bool ASMs3DLag::integrateEdge (Integrand& integrand, int lEdge,
if (skipMe) continue;
if (lEdge < 5)
ip = i1*nGauss;
ip = i1*ng;
else if (lEdge < 9)
ip = i2*nGauss;
ip = i2*ng;
else
ip = i3*nGauss;
ip = i3*ng;
// Set up nodal point coordinates for current element
if (!this->getElementCoordinates(Xnod,iel)) return false;
@ -792,7 +802,7 @@ bool ASMs3DLag::integrateEdge (Integrand& integrand, int lEdge,
fe.iGP = firstp + ip; // Global integration point counter
for (int i = 0; i < nGauss && ok; i++, fe.iGP++)
for (int i = 0; i < ng && ok; i++, fe.iGP++)
{
// Gauss point coordinates on the edge
xi[lDir] = xg[i];
@ -836,14 +846,13 @@ int ASMs3DLag::evalPoint (const double* xi, double* param, Vec3& X) const
{
// Evaluate the parametric values of the point and nodes
std::array<RealArray,3> u;
std::array<int,3> p({{p1,p2,p3}});
for (int d = 0; d < 3; d++)
{
if (svol)
param[d] = (1.0-xi[d])*svol->startparam(d) + xi[d]*svol->endparam(d);
else
param[d] = xi[d];
if (!this->getGridParameters(u[d],d,p[d]-1)) return -3;
if (!this->getGridParameters(u[d],d,svol->order(d)-1)) return -3;
}
// Search for the closest node

25
src/ASM/ASMs3Drecovery.C
View File

@ -175,14 +175,17 @@ bool ASMs3D::assembleL2matrices (SparseMatrix& A, StdVector& B,
const int nel2 = n2 - p2 + 1;
const int nel3 = n3 - p3 + 1;
int pmax = p1 > p2 ? p1 : p2;
if (pmax < p3) pmax = p3;
// Get Gaussian quadrature point coordinates (and weights if continuous)
const int ng1 = continuous ? nGauss : p1 - 1;
const int ng2 = continuous ? nGauss : p2 - 1;
const int ng3 = continuous ? nGauss : p3 - 1;
const int ng1 = continuous ? this->getNoGaussPt(pmax,true) : p1 - 1;
const int ng2 = continuous ? ng1 : p2 - 1;
const int ng3 = continuous ? ng2 : p3 - 1;
const double* xg = GaussQuadrature::getCoord(ng1);
const double* yg = GaussQuadrature::getCoord(ng2);
const double* zg = GaussQuadrature::getCoord(ng3);
const double* wg = continuous ? GaussQuadrature::getWeight(nGauss) : 0;
const double* wg = continuous ? GaussQuadrature::getWeight(ng1) : 0;
if (!xg || !yg || !zg) return false;
if (continuous && !wg) return false;
@ -346,26 +349,30 @@ Go::SplineVolume* ASMs3D::projectSolutionLeastSquare (const IntegrandBase& integ
PROFILE1("test L2- projection");
// Compute parameter values of the result sampling points (Gauss-Interpl. points)
// Get Gaussian quadrature points and weights
const double* xg = GaussQuadrature::getCoord(nGauss);
const double* wg = GaussQuadrature::getWeight(nGauss);
int p = svol->order(0);
for (int d = 1; d < 3; d++)
if (p < svol->order(d)) p = svol->order(d);
const int ng = this->getNoGaussPt(p,true);
const double* xg = GaussQuadrature::getCoord(ng);
const double* wg = GaussQuadrature::getWeight(ng);
if (!xg || !wg) return nullptr;
std::array<Matrix,3> ggpar;
std::array<RealArray,3> gpar, wgpar;
for (int dir = 0; dir < 3; dir++)
{
this->getGaussPointParameters(ggpar[dir],dir,nGauss,xg);
this->getGaussPointParameters(ggpar[dir],dir,ng,xg);
gpar[dir] = ggpar[dir];
// Gauss weights at parameter values
const Go::BsplineBasis& basis = svol->basis(dir);
RealArray::const_iterator knotit = basis.begin();
RealArray& tmp = wgpar[dir];
tmp.reserve(nGauss*(basis.numCoefs()-basis.order()));
tmp.reserve(ng*(basis.numCoefs()-basis.order()));
for (int i = basis.order(); i <= basis.numCoefs(); i++)
{
double d = knotit[i]-knotit[i-1];
for (int j = 0; j < nGauss; j++)
for (int j = 0; j < ng; j++)
tmp.push_back(d > 0.0 ? wg[j]*d*0.5 : 0.0);
}
}

3
src/ASM/FiniteElement.C
View File

@ -12,6 +12,7 @@
//==============================================================================
#include "FiniteElement.h"
#include "Vec3Oper.h"
std::ostream& FiniteElement::write (std::ostream& os) const
@ -21,6 +22,8 @@ std::ostream& FiniteElement::write (std::ostream& os) const
<<"\n u, v, w: "<< u <<" "<< v <<" "<< w
<<"\n xi, eta, zeta: "<< xi <<" "<< eta <<" "<< zeta
<<"\n h, detJxW: "<< h <<" "<< detJxW << std::endl;
for (size_t n = 0; n < XC.size(); n++)
os <<"\n XC_"<< n+1 <<": "<< XC[n];
if (!N.empty()) os <<"N:"<< N;
if (!dNdX.empty()) os <<"dNdX:"<< dNdX;
if (!d2NdX2.empty()) os <<"d2NdX2:"<< d2NdX2;