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

Fixed evaluation of secondary solution when integrandType==2: Must use 2nd derivatives

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1202 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo 2011-09-24 17:53:03 +00:00 committed by Knut Morten Okstad
parent 56e2e2642f
commit e52597abb6
2 changed files with 82 additions and 102 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -692,7 +692,7 @@ Vec3 ASMu2D::getCoord (size_t inod) const
} }
void ASMu2D::getGaussPointParameters (Vector& uGP, int dir, int nGauss, void ASMu2D::getGaussPointParameters (RealArray& uGP, int dir, int nGauss,
int iel, const double* xi) const int iel, const double* xi) const
{ {
#ifdef INDEX_CHECK #ifdef INDEX_CHECK
@ -709,8 +709,8 @@ void ASMu2D::getGaussPointParameters (Vector& uGP, int dir, int nGauss,
double ustart = (dir==0) ? el->umin() : el->vmin(); double ustart = (dir==0) ? el->umin() : el->vmin();
double ustop = (dir==0) ? el->umax() : el->vmax(); double ustop = (dir==0) ? el->umax() : el->vmax();
for (int i = 1; i <= nGauss; i++) for (int i = 0; i < nGauss; i++)
uGP(i) = 0.5*((ustop-ustart)*xi[i-1] + ustop+ustart); uGP[i] = 0.5*((ustop-ustart)*xi[i] + ustop+ustart);
} }
@ -836,7 +836,7 @@ bool ASMu2D::integrate (Integrand& integrand,
if (!this->getElementCoordinates(Xnod,iel)) return false; if (!this->getElementCoordinates(Xnod,iel)) return false;
// Compute parameter values of the Gauss points over this element // Compute parameter values of the Gauss points over this element
Vector gpar[2], redpar[2]; RealArray gpar[2], redpar[2];
for (int d = 0; d < 2; d++) for (int d = 0; d < 2; d++)
{ {
this->getGaussPointParameters(gpar[d],d,nGauss,iel,xg); this->getGaussPointParameters(gpar[d],d,nGauss,iel,xg);
@ -899,25 +899,25 @@ bool ASMu2D::integrate (Integrand& integrand,
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1]; LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
#if 0
if (integrand.getIntegrandType() > 10) if (integrand.getIntegrandType() > 10)
{ {
// --- Selective reduced integration loop ------------------------------ // --- Selective reduced integration loop ------------------------------
int ip = ((i2-p2)*nRed*nel1 + i1-p1)*nRed; for (int j = 0; j < nRed; j++)
for (int j = 0; j < nRed; j++, ip += nRed*(nel1-1)) for (int i = 0; i < nRed; i++)
for (int i = 0; i < nRed; i++, ip++)
{ {
// Local element coordinates of current integration point // Local element coordinates of current integration point
fe.xi = xr[i]; fe.xi = xr[i];
fe.eta = xr[j]; fe.eta = xr[j];
// Parameter values of current integration point // Parameter values of current integration point
fe.u = redpar[0](i+1,i1-p1+1); fe.u = redpar[0][i];
fe.v = redpar[1](j+1,i2-p2+1); fe.v = redpar[1][j];
// Fetch basis function derivatives at current point // Compute basis function derivatives at current point
extractBasis(splineRed[ip],fe.N,dNdu); Go::BasisDerivsSf spline;
lrspline->computeBasis(fe.u,fe.v,spline,iel-1);
extractBasis(spline,fe.N,dNdu);
// Compute Jacobian inverse and derivatives // Compute Jacobian inverse and derivatives
fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu);
@ -927,7 +927,6 @@ bool ASMu2D::integrate (Integrand& integrand,
return false; return false;
} }
} }
#endif
// --- Integration loop over all Gauss points in each direction ---------- // --- Integration loop over all Gauss points in each direction ----------
@ -940,35 +939,28 @@ bool ASMu2D::integrate (Integrand& integrand,
fe.eta = xg[j]; fe.eta = xg[j];
// Parameter values of current integration point // Parameter values of current integration point
fe.u = gpar[0](i+1); fe.u = gpar[0][i];
fe.v = gpar[1](j+1); fe.v = gpar[1][j];
// compute basis function values at integration point // Compute basis function derivatives at current integration point
Go::BasisDerivsSf spline; if (integrand.getIntegrandType() == 2) {
Go::BasisDerivsSf2 spline2; Go::BasisDerivsSf2 spline;
Go::BasisDerivsSf splineRed; lrspline->computeBasis(fe.u,fe.v,spline, iel-1);
if (integrand.getIntegrandType() == 2) extractBasis(spline,fe.N,dNdu,d2Ndu2);
lrspline->computeBasis(fe.u,fe.v,spline2, iel-1); }
else else {
Go::BasisDerivsSf spline;
lrspline->computeBasis(fe.u,fe.v,spline, iel-1); lrspline->computeBasis(fe.u,fe.v,spline, iel-1);
if (integrand.getIntegrandType() > 10)
lrspline->computeBasis(fe.u,fe.v,splineRed, iel-1);
// Fetch basis function derivatives at current integration point
if (integrand.getIntegrandType() == 2)
extractBasis(spline2,fe.N,dNdu,d2Ndu2);
else
extractBasis(spline,fe.N,dNdu); extractBasis(spline,fe.N,dNdu);
#if SP_DEBUG > 4 #if SP_DEBUG > 4
std::cout <<"\nBasis functions at a integration point "; std::cout <<"\nBasis functions at a integration point ";
std::cout <<" : (u,v) = "<< spline.param[0] <<" "<< spline.param[1] std::cout <<" : (u,v) = "<< spline.param[0] <<" "<< spline.param[1]
<<" left_idx = "<< spline.left_idx[0] <<" "<< spline.left_idx[1] << std::endl; <<" left_idx = "<< spline.left_idx[0] <<" "<< spline.left_idx[1] << std::endl;
for (size_t ii = 0; ii < spline.basisValues.size(); ii++) for (size_t ii = 0; ii < spline.basisValues.size(); ii++)
std::cout << 1+ii <<'\t'<< spline.basisValues[ii] <<'\t' std::cout << 1+ii <<'\t'<< spline.basisValues[ii] <<'\t'
<< spline.basisDerivs_u[ii] <<'\t'<< spline.basisDerivs_v[ii] << std::endl; << spline.basisDerivs_u[ii] <<'\t'<< spline.basisDerivs_v[ii] << std::endl;
#endif #endif
}
// Compute Jacobian inverse of coordinate mapping and derivatives // Compute Jacobian inverse of coordinate mapping and derivatives
fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu); fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu);
@ -1277,21 +1269,6 @@ bool ASMu2D::tesselate (ElementBlock& grid, const int* npe) const
return true; return true;
} }
#if 0
void ASMu2D::scatterInd (int n1, int n2, int p1, int p2,
const int* start, IntVec& index)
{
index.reserve(p1*p2);
int ip = ((start[1]-p2+1))*n1 + (start[0]-p1+1);
for (int i2 = 0; i2 < p2; i2++, ip += n1-p1)
for (int i1 = 0; i1 < p1; i1++, ip++)
index.push_back(ip);
}
#endif
bool ASMu2D::evalSolution (Matrix& sField, const Vector& locSol, bool ASMu2D::evalSolution (Matrix& sField, const Vector& locSol,
const int* npe) const const int* npe) const
@ -1335,7 +1312,7 @@ bool ASMu2D::evalSolution (Matrix& sField, const Vector& locSol,
int iel = i/4; // points are always listed in the same order as the elemnts, 4 pts per element int iel = i/4; // points are always listed in the same order as the elemnts, 4 pts per element
// fetch index of non-zero basis functions on this element // fetch index of non-zero basis functions on this element
IntVec ip = MNPC[iel]; const IntVec& ip = MNPC[iel];
// evaluate the basis functions at current parametric point // evaluate the basis functions at current parametric point
Go::BasisPtsSf spline; Go::BasisPtsSf spline;
@ -1434,56 +1411,65 @@ bool ASMu2D::evalSolution (Matrix& sField, const Integrand& integrand,
const RealArray* gpar, bool regular) const const RealArray* gpar, bool regular) const
{ {
#ifdef SP_DEBUG #ifdef SP_DEBUG
std::cout << "ASMu2D::evalSolution( )\n"; std::cout <<"ASMu2D::evalSolution(Matrix&,const Integrand&,const RealArray*,bool)\n";
#endif #endif
// TODO: investigate the possibility of doing "regular" refinement by unfiorm
// tesselation grid and ignoring LR meshlines
if(gpar[0].size() != gpar[1].size()) sField.resize(0,0);
return false;
sField.resize(0,0); // TODO: investigate the possibility of doing "regular" refinement by
// uniform tesselation grid and ignoring LR mesh lines
// Fetch nodal (control point) coordinates if (gpar[0].size() != gpar[1].size())
Matrix Xnod; return false;
IntVec ip; Vector solPt;
Vector solPt; Matrix dNdu, dNdX, Jac, Xnod;
Matrix dNdu, dNdX, Jac; Matrix3D d2Ndu2, d2NdX2, Hess;
// Evaluate the secondary solution field at each point size_t nPoints = gpar[0].size();
size_t nPoints = gpar[0].size(); bool use2ndDer = integrand.getIntegrandType() == 2;
for (size_t i = 0; i < nPoints; i++) // Evaluate the secondary solution field at each point
{ for (size_t i = 0; i < nPoints; i++)
// fetch element containing evaluation point {
int iel = i/4; // points are always listed in the same order as the elemnts, 4 pts per element // Fetch element containing evaluation point. Points are always listed
int iel = i/4; // in the same order as the elements, 4 pts per element.
// fetch number of nonzero basis functions on this element as well as their indices // Evaluate the basis functions at current parametric point
int nBasis = lrspline->getElement(iel)->nBasisFunctions(); Vector N(lrspline->getElement(iel)->nBasisFunctions());
Vector N(nBasis); if (use2ndDer)
ip = MNPC[iel]; {
Go::BasisDerivsSf2 spline;
lrspline->computeBasis(gpar[0][i],gpar[1][i],spline,iel);
extractBasis(spline,N,dNdu,d2Ndu2);
}
else
{
Go::BasisDerivsSf spline;
lrspline->computeBasis(gpar[0][i],gpar[1][i],spline,iel);
extractBasis(spline,N,dNdu);
}
// evaluate the basis functions at current parametric point // Set up control point (nodal) coordinates for current element
Go::BasisDerivsSf spline; if (!this->getElementCoordinates(Xnod,iel+1)) return false;
lrspline->computeBasis(gpar[0][i], gpar[1][i], spline, iel);
extractBasis(spline, N, dNdu);
// Set up control point (nodal) coordinates for current element // Compute the Jacobian inverse
if (!this->getElementCoordinates(Xnod,iel+1)) return false; if (utl::Jacobian(Jac,dNdX,Xnod,dNdu) == 0.0) // Jac = (Xnod * dNdu)^-1
continue; // skip singular points
// Compute the Jacobian inverse // Compute Hessian of coordinate mapping and 2nd order derivatives
if (utl::Jacobian(Jac,dNdX,Xnod,dNdu) == 0.0) // Jac = (Xnod * dNdu)^-1 if (use2ndDer)
continue; // skip singular points if (!utl::Hessian(Hess,d2NdX2,Jac,Xnod,d2Ndu2,dNdu))
continue;
// Now evaluate the solution field // Now evaluate the solution field
if (!integrand.evalSol(solPt,N,dNdX,Xnod*N,ip)) if (!integrand.evalSol(solPt,N,dNdX,d2NdX2,Xnod*N,MNPC[iel]))
return false; return false;
else if (sField.empty()) else if (sField.empty())
sField.resize(solPt.size(),nPoints,true); sField.resize(solPt.size(),nPoints,true);
sField.fillColumn(1+i,solPt); sField.fillColumn(1+i,solPt);
} }
return true; return true;
} }

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

@ -149,13 +149,13 @@ public:
//! \param[in] dir Parameter direction to refine //! \param[in] dir Parameter direction to refine
//! \param[in] nInsert Number of extra knots to insert in each knot-span //! \param[in] nInsert Number of extra knots to insert in each knot-span
bool uniformRefine(int dir, int nInsert); bool uniformRefine(int dir, int nInsert);
//! \brief Refine the parametrization by inserting extra tensor knots. //! \brief Refines the parametrization by inserting extra tensor knots.
//! \details This method is mainly kept for backward compatability with the //! \details This method is mainly kept for backward compatability with the
//! "REFINE" keyword in the input file. //! "REFINE" keyword in the input file.
//! \param[in] dir Parameter direction to refine //! \param[in] dir Parameter direction to refine
//! \param[in] xi Relative positions of added knots in each existing knot span //! \param[in] xi Relative positions of added knots in each existing knot span
bool refine(int dir, const RealArray& xi); bool refine(int dir, const RealArray& xi);
//! \brief Raise the order of the tensor spline object for this patch. //! \brief Raises the order of the tensor spline object for this patch.
//! \details This method is mainly kept for backward compatability with the //! \details This method is mainly kept for backward compatability with the
//! "RAISEORDER" keyword in the input file. //! "RAISEORDER" keyword in the input file.
//! \param[in] ru Number of times to raise the order in u-direction //! \param[in] ru Number of times to raise the order in u-direction
@ -352,7 +352,7 @@ protected:
//! \param[in] nGauss Number of Gauss points along a knot-span //! \param[in] nGauss Number of Gauss points along a knot-span
//! \param[in] iel Element index //! \param[in] iel Element index
//! \param[in] xi Dimensionless Gauss point coordinates [-1,1] //! \param[in] xi Dimensionless Gauss point coordinates [-1,1]
void getGaussPointParameters(Vector& uGP, int dir, int nGauss, void getGaussPointParameters(RealArray& uGP, int dir, int nGauss,
int iel, const double* xi) const; int iel, const double* xi) const;
//! \brief Calculates parameter values for the Greville points. //! \brief Calculates parameter values for the Greville points.
@ -384,16 +384,10 @@ protected:
static void extractBasis(const Go::BasisDerivsSf2& spline, static void extractBasis(const Go::BasisDerivsSf2& spline,
Vector& N, Matrix& dNdu, Matrix3D& d2Ndu2); Vector& N, Matrix& dNdu, Matrix3D& d2Ndu2);
//! \brief Auxilliary function for computation of basis function indices.
static void scatterInd(int n1, int n2, int p1, int p2,
const int* start, IntVec& index);
private:
protected: protected:
LR::LRSplineSurface* lrspline; //!< Pointer to the actual spline surface object LR::LRSplineSurface* lrspline; //!< Pointer to the actual spline surface object
Go::SplineSurface* tensorspline; //!< Pointer to the original tensor spline object private:
Go::SplineSurface* tensorspline; //!< Pointer to the original tensor spline object
// the tensor spline object is kept for backwards compatability with the REFINE and RAISEORDER // the tensor spline object is kept for backwards compatability with the REFINE and RAISEORDER
// key-words, although we take note that there is a possibility of optimization since all mapping // key-words, although we take note that there is a possibility of optimization since all mapping
// values and jacobians may be performed on this object for increased efficiency // values and jacobians may be performed on this object for increased efficiency