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

Added: Support for third order derivatives in ASMu2D

Browse Source
This commit is contained in:
Arne Morten Kvarving 2018-11-24 23:14:35 +01:00 committed by Knut Morten Okstad
parent b6b6c2dd7d
commit 712a848874
1 changed files with 32 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -968,6 +968,9 @@ bool ASMu2D::integrate (Integrand& integrand,
PROFILE2("ASMu2D::integrate(I)");
bool use2ndDer = integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES;
bool use3rdDer = integrand.getIntegrandType() & Integrand::THIRD_DERIVATIVES;
// Get Gaussian quadrature points and weights
const double* xg = GaussQuadrature::getCoord(nGauss);
const double* wg = GaussQuadrature::getWeight(nGauss);
@ -992,8 +995,11 @@ bool ASMu2D::integrate (Integrand& integrand,
std::vector<Go::BasisDerivsSf> spline1, splineRed;
std::vector<Go::BasisDerivsSf2> spline2;
std::vector<Go::BasisDerivsSf3> spline3;
if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES)
if (use3rdDer)
spline3.resize(nel*nGauss*nGauss);
else if (use2ndDer)
spline2.resize(nel*nGauss*nGauss);
else
spline1.resize(nel*nGauss*nGauss);
@ -1008,7 +1014,9 @@ bool ASMu2D::integrate (Integrand& integrand,
this->getGaussPointParameters(v,1,nGauss,1+iel,xg);
for (int j = 0; j < nGauss; j++)
for (int i = 0; i < nGauss; i++, jp++)
if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES)
if (use3rdDer)
lrspline->computeBasis(u[i],v[j],spline3[jp],iel);
else if (use2ndDer)
lrspline->computeBasis(u[i],v[j],spline2[jp],iel);
else
lrspline->computeBasis(u[i],v[j],spline1[jp],iel);
@ -1050,6 +1058,7 @@ bool ASMu2D::integrate (Integrand& integrand,
fe.q = lrspline->order(1) - 1;
Matrix dNdu, Xnod, Jac;
Matrix3D d2Ndu2, Hess;
Matrix4D d3Ndu3;
double dXidu[2];
double param[3] = { 0.0, 0.0, 0.0 };
Vec4 X(param);
@ -1173,7 +1182,9 @@ bool ASMu2D::integrate (Integrand& integrand,
fe.v = param[1] = gpar[1][j];
// Extract basis function derivatives at current integration point
if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES)
if (use3rdDer)
SplineUtils::extractBasis(spline3[fe.iGP-firstIp],fe.N,dNdu,d2Ndu2,d3Ndu3);
else if (use2ndDer)
SplineUtils::extractBasis(spline2[fe.iGP-firstIp],fe.N,dNdu,d2Ndu2);
else {
SplineUtils::extractBasis(spline1[fe.iGP-firstIp],fe.N,dNdu);
@ -1199,7 +1210,7 @@ bool ASMu2D::integrate (Integrand& integrand,
if (fe.detJxW == 0.0) continue; // skip singular points
// Compute Hessian of coordinate mapping and 2nd order derivatives
if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES)
if (use2ndDer)
{
if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu))
ok = false;
@ -1207,6 +1218,10 @@ bool ASMu2D::integrate (Integrand& integrand,
utl::Hessian(Hess,fe.H);
}
// Compute 3rd order derivatives
if (use3rdDer)
ok &= utl::Hessian2(fe.d3NdX3,Jac,d3Ndu3);
// Compute G-matrix
if (integrand.getIntegrandType() & Integrand::G_MATRIX)
utl::getGmat(Jac,dXidu,fe.G);
@ -1887,6 +1902,7 @@ bool ASMu2D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
size_t nPoints = gpar[0].size();
bool use2ndDer = integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES;
bool use3rdDer = integrand.getIntegrandType() & Integrand::THIRD_DERIVATIVES;
if (nPoints != gpar[1].size())
return false;
@ -1896,6 +1912,7 @@ bool ASMu2D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
Vector solPt;
Matrix dNdu, Jac, Xnod;
Matrix3D d2Ndu2, Hess;
Matrix4D d3Ndu3;
// Evaluate the secondary solution field at each point
for (size_t i = 0; i < nPoints; i++, fe.iGP++)
@ -1907,7 +1924,13 @@ bool ASMu2D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
int iel = lrspline->getElementContaining(fe.u,fe.v);
// Evaluate the basis functions at current parametric point
if (use2ndDer)
if (use3rdDer)
{
Go::BasisDerivsSf3 spline;
lrspline->computeBasis(gpar[0][i],gpar[1][i],spline,iel);
SplineUtils::extractBasis(spline,fe.N,dNdu,d2Ndu2,d3Ndu3);
}
else if (use2ndDer)
{
Go::BasisDerivsSf2 spline;
lrspline->computeBasis(fe.u,fe.v,spline,iel);
@ -1935,6 +1958,10 @@ bool ASMu2D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
utl::Hessian(Hess,fe.H);
}
// Compute 3rd order derivatives
if (use3rdDer)
utl::Hessian2(fe.d3NdX3,Jac,d3Ndu3);
// Store tangent vectors in fe.G for shells
if (nsd > 2) fe.G = Jac;