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

Fixed: Account for that the number of space dimensions in a patch

Browse Source 
can be less than the number of coordinates on file
This commit is contained in:
Knut Morten Okstad 2018-02-10 17:16:22 +01:00
parent f9754b1912
commit ac6c1c3a6b
13 changed files with 102 additions and 212 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/ASM/ASMs1D.C
View File

@ -686,7 +686,7 @@ double ASMs1D::getElementEnds (int i, Vec3Vec& XC) const
XC.reserve(elmCS.empty() ? 2 : 3); XC.reserve(elmCS.empty() ? 2 : 3);
const double* pt = &XYZ.front(); const double* pt = &XYZ.front();
for (int j = 0; j < 2; j++, pt += dim) for (int j = 0; j < 2; j++, pt += dim)
XC.push_back(Vec3(pt,dim)); XC.push_back(Vec3(pt,nsd));
// Calculate the characteristic element size // Calculate the characteristic element size
double h = getElementSize(XC); double h = getElementSize(XC);

2
src/ASM/ASMs2D.C
View File

@ -1456,7 +1456,7 @@ double ASMs2D::getElementCorners (int i1, int i2, Vec3Vec& XC) const
XC.reserve(4); XC.reserve(4);
const double* pt = &XYZ.front(); const double* pt = &XYZ.front();
for (int i = 0; i < 4; i++, pt += dim) for (int i = 0; i < 4; i++, pt += dim)
XC.push_back(Vec3(pt,dim)); XC.push_back(Vec3(pt,nsd));
return getElementSize(XC); return getElementSize(XC);
} }

2
src/ASM/ASMs3D.C
View File

@ -1659,7 +1659,7 @@ double ASMs3D::getElementCorners (int i1, int i2, int i3, Vec3Vec& XC) const
XC.reserve(8); XC.reserve(8);
const double* pt = &XYZ.front(); const double* pt = &XYZ.front();
for (int i = 0; i < 8; i++, pt += dim) for (int i = 0; i < 8; i++, pt += dim)
XC.push_back(Vec3(pt,dim)); XC.push_back(Vec3(pt,nsd));
return getElementSize(XC); return getElementSize(XC);
} }

51
src/ASM/SplineField2D.C
View File

@ -35,6 +35,8 @@ SplineField2D::SplineField2D (const ASMs2D* patch,
const int p2 = basis->order_v(); const int p2 = basis->order_v();
nelm = (n1-p1+1)*(n2-p2+1); nelm = (n1-p1+1)*(n2-p2+1);
nsd = patch->getNoSpaceDim();
// Ensure the values array has compatible length, pad with zeros if necessary // Ensure the values array has compatible length, pad with zeros if necessary
size_t ofs = 0; size_t ofs = 0;
for (char i = 1; i < nbasis; ++i) for (char i = 1; i < nbasis; ++i)
@ -93,7 +95,7 @@ double SplineField2D::valueCoor (const Vec4& x) const
pt[1] = x[1]; pt[1] = x[1];
pt[2] = x[2]; pt[2] = x[2];
double clo_u, clo_v, dist; double clo_u, clo_v, dist;
#pragma omp critical #pragma omp critical
surf->closestPoint(pt, clo_u, clo_v, clopt, dist, 1e-10); surf->closestPoint(pt, clo_u, clo_v, clopt, dist, 1e-10);
fe.u = clo_u; fe.u = clo_u;
fe.v = clo_v; fe.v = clo_v;
@ -184,9 +186,9 @@ bool SplineField2D::gradFE (const FiniteElement& fe, Vector& grad) const
uorder,vorder,spline.left_idx,ip); uorder,vorder,spline.left_idx,ip);
// Evaluate the Jacobian inverse // Evaluate the Jacobian inverse
Matrix Xnod, Jac; Matrix Xnod(nsd,ip.size()), Jac;
Vector Xctrl(&(*surf->coefs_begin()),surf->coefs_end()-surf->coefs_begin()); for (size_t i = 0; i < ip.size(); i++)
utl::gather(ip,surf->dimension(),Xctrl,Xnod); Xnod.fillColumn(1+i,&(*surf->coefs_begin())+surf->dimension()*ip[i]);
utl::Jacobian(Jac,dNdX,Xnod,dNdu); utl::Jacobian(Jac,dNdX,Xnod,dNdu);
// Evaluate the gradient of the solution field at the given point // Evaluate the gradient of the solution field at the given point
@ -222,73 +224,38 @@ bool SplineField2D::hessianFE (const FiniteElement& fe, Matrix& H) const
if (!basis) return false; if (!basis) return false;
if (!surf) return false; if (!surf) return false;
// Order of basis
const int uorder = surf->order_u();
const int vorder = surf->order_v();
const size_t nen = uorder*vorder;
// Evaluate the basis functions at the given point
Go::BasisDerivsSf spline;
Go::BasisDerivsSf2 spline2; Go::BasisDerivsSf2 spline2;
Matrix3D d2Ndu2; Matrix3D d2Ndu2;
Matrix dNdu, dNdX;
IntVec ip; IntVec ip;
if (surf == basis) { if (surf == basis) {
#pragma omp critical #pragma omp critical
surf->computeBasis(fe.u,fe.v,spline2); surf->computeBasis(fe.u,fe.v,spline2);
dNdu.resize(nen,2); const size_t nen = surf->order_u()*surf->order_v();
d2Ndu2.resize(nen,2,2); d2Ndu2.resize(nen,2,2);
for (size_t n = 1; n <= nen; n++) { for (size_t n = 1; n <= nen; n++) {
dNdu(n,1) = spline2.basisDerivs_u[n-1];
dNdu(n,2) = spline2.basisDerivs_v[n-1];
d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1]; d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1];
d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1]; d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1];
d2Ndu2(n,2,2) = spline2.basisDerivs_vv[n-1]; d2Ndu2(n,2,2) = spline2.basisDerivs_vv[n-1];
} }
ASMs2D::scatterInd(surf->numCoefs_u(),surf->numCoefs_v(), ASMs2D::scatterInd(surf->numCoefs_u(),surf->numCoefs_v(),
uorder,vorder,spline2.left_idx,ip); surf->order_u(),surf->order_v(),
spline2.left_idx,ip);
} }
else { else {
#pragma omp critical
surf->computeBasis(fe.u,fe.v,spline);
dNdu.resize(nen,2);
for (size_t n = 1; n <= nen; n++) {
dNdu(n,1) = spline.basisDerivs_u[n-1];
dNdu(n,2) = spline.basisDerivs_v[n-1];
}
ASMs2D::scatterInd(surf->numCoefs_u(),surf->numCoefs_v(),
uorder,vorder,spline.left_idx,ip);
}
// Evaluate the Jacobian inverse
Matrix Xnod, Jac;
Vector Xctrl(&(*surf->coefs_begin()),surf->coefs_end()-surf->coefs_begin());
utl::gather(ip,surf->dimension(),Xctrl,Xnod);
utl::Jacobian(Jac,dNdX,Xnod,dNdu);
// Evaluate the gradient of the solution field at the given point
if (basis != surf)
{
// Mixed formulation, the solution uses a different basis than the geometry // Mixed formulation, the solution uses a different basis than the geometry
#pragma omp critical #pragma omp critical
basis->computeBasis(fe.u,fe.v,spline2); basis->computeBasis(fe.u,fe.v,spline2);
const size_t nbf = basis->order_u()*basis->order_v(); const size_t nbf = basis->order_u()*basis->order_v();
dNdu.resize(nbf,2);
d2Ndu2.resize(nbf,2,2); d2Ndu2.resize(nbf,2,2);
for (size_t n = 1; n <= nbf; n++) { for (size_t n = 1; n <= nbf; n++) {
dNdu(n,1) = spline2.basisDerivs_u[n-1];
dNdu(n,2) = spline2.basisDerivs_v[n-1];
d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1]; d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1];
d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1]; d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1];
d2Ndu2(n,2,2) = spline2.basisDerivs_vv[n-1]; d2Ndu2(n,2,2) = spline2.basisDerivs_vv[n-1];
} }
ip.clear();
ASMs2D::scatterInd(basis->numCoefs_u(),basis->numCoefs_v(), ASMs2D::scatterInd(basis->numCoefs_u(),basis->numCoefs_v(),
basis->order_u(),basis->order_v(), basis->order_u(),basis->order_v(),
spline2.left_idx,ip); spline2.left_idx,ip);

3
src/ASM/SplineField2D.h
View File

@ -77,6 +77,9 @@ public:
protected: protected:
const Go::SplineSurface* basis; //!< Spline basis description const Go::SplineSurface* basis; //!< Spline basis description
const Go::SplineSurface* surf; //!< Spline geometry description const Go::SplineSurface* surf; //!< Spline geometry description
private:
unsigned char nsd; //!< Number of space dimensions
}; };
#endif #endif

53
src/ASM/SplineField3D.C
View File

@ -37,6 +37,8 @@ SplineField3D::SplineField3D (const ASMs3D* patch,
const int p3 = basis->order(2); const int p3 = basis->order(2);
nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1); nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1);
nsd = patch->getNoSpaceDim();
size_t ofs = 0; size_t ofs = 0;
for (char i = 1; i < nbasis; ++i) for (char i = 1; i < nbasis; ++i)
ofs += patch->getNoNodes(i)*patch->getNoFields(i); ofs += patch->getNoNodes(i)*patch->getNoFields(i);
@ -95,7 +97,7 @@ double SplineField3D::valueCoor (const Vec4& x) const
pt[1] = x[1]; pt[1] = x[1];
pt[2] = x[2]; pt[2] = x[2];
double clo_u, clo_v, clo_w, dist; double clo_u, clo_v, clo_w, dist;
#pragma omp critical #pragma omp critical
vol->closestPoint(pt, clo_u, clo_v, clo_w, clopt, dist, 1e-5); vol->closestPoint(pt, clo_u, clo_v, clo_w, clopt, dist, 1e-5);
fe.u = clo_u; fe.u = clo_u;
@ -191,9 +193,9 @@ bool SplineField3D::gradFE (const FiniteElement& fe, Vector& grad) const
uorder,vorder,worder,spline.left_idx,ip); uorder,vorder,worder,spline.left_idx,ip);
// Evaluate the Jacobian inverse // Evaluate the Jacobian inverse
Matrix Xnod, Jac; Matrix Xnod(nsd,ip.size()), Jac;
Vector Xctrl(&(*vol->coefs_begin()),vol->coefs_end()-vol->coefs_begin()); for (size_t i = 0; i < ip.size(); i++)
utl::gather(ip,vol->dimension(),Xctrl,Xnod); Xnod.fillColumn(1+i,&(*vol->coefs_begin())+vol->dimension()*ip[i]);
utl::Jacobian(Jac,dNdX,Xnod,dNdu); utl::Jacobian(Jac,dNdX,Xnod,dNdu);
// Evaluate the gradient of the solution field at the given point // Evaluate the gradient of the solution field at the given point
@ -230,25 +232,16 @@ bool SplineField3D::hessianFE (const FiniteElement& fe, Matrix& H) const
if (!basis) return false; if (!basis) return false;
if (!vol) return false; if (!vol) return false;
const int uorder = vol->order(0);
const int vorder = vol->order(1);
const int worder = vol->order(2);
const size_t nen = uorder*vorder*worder;
// Evaluate the basis functions at the given point
Go::BasisDerivs spline;
Go::BasisDerivs2 spline2; Go::BasisDerivs2 spline2;
Matrix3D d2Ndu2; Matrix3D d2Ndu2;
Matrix dNdu(nen,3), dNdX;
IntVec ip; IntVec ip;
if (vol == basis) { if (vol == basis) {
#pragma omp critical #pragma omp critical
vol->computeBasis(fe.u,fe.v,fe.w,spline2); vol->computeBasis(fe.u,fe.v,fe.w,spline2);
const size_t nen = vol->order(0)*vol->order(1)*vol->order(2);
d2Ndu2.resize(nen,3,3); d2Ndu2.resize(nen,3,3);
for (size_t n = 1; n <= nen; n++) { for (size_t n = 1; n <= nen; n++) {
dNdu(n,1) = spline2.basisDerivs_u[n-1];
dNdu(n,2) = spline2.basisDerivs_v[n-1];
dNdu(n,3) = spline2.basisDerivs_w[n-1];
d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1]; d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1];
d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1]; d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1];
d2Ndu2(n,1,3) = d2Ndu2(n,3,1) = spline2.basisDerivs_uw[n-1]; d2Ndu2(n,1,3) = d2Ndu2(n,3,1) = spline2.basisDerivs_uw[n-1];
@ -256,40 +249,19 @@ bool SplineField3D::hessianFE (const FiniteElement& fe, Matrix& H) const
d2Ndu2(n,2,3) = d2Ndu2(n,3,2) = spline2.basisDerivs_vw[n-1]; d2Ndu2(n,2,3) = d2Ndu2(n,3,2) = spline2.basisDerivs_vw[n-1];
d2Ndu2(n,3,3) = spline2.basisDerivs_ww[n-1]; d2Ndu2(n,3,3) = spline2.basisDerivs_ww[n-1];
} }
ASMs3D::scatterInd(vol->numCoefs(0),vol->numCoefs(1),vol->numCoefs(2), ASMs3D::scatterInd(vol->numCoefs(0),vol->numCoefs(1),vol->numCoefs(2),
uorder,vorder,worder,spline2.left_idx,ip); vol->order(0),vol->order(1),vol->order(2),
spline2.left_idx,ip);
} }
else { else {
#pragma omp critical
vol->computeBasis(fe.u,fe.v,fe.w,spline);
for (size_t n = 1; n <= nen; n++) {
dNdu(n,1) = spline.basisDerivs_u[n-1];
dNdu(n,2) = spline.basisDerivs_v[n-1];
dNdu(n,3) = spline.basisDerivs_w[n-1];
}
ASMs3D::scatterInd(vol->numCoefs(0),vol->numCoefs(1),vol->numCoefs(2),
uorder,vorder,worder,spline.left_idx,ip);
}
// Evaluate the Jacobian inverse
Matrix Xnod, Jac;
Vector Xctrl(&(*vol->coefs_begin()),vol->coefs_end()-vol->coefs_begin());
utl::gather(ip,vol->dimension(),Xctrl,Xnod);
utl::Jacobian(Jac,dNdX,Xnod,dNdu);
// Evaluate the gradient of the solution field at the given point
if (basis != vol) {
// Mixed formulation, the solution uses a different basis than the geometry // Mixed formulation, the solution uses a different basis than the geometry
#pragma omp critical #pragma omp critical
basis->computeBasis(fe.u,fe.v,fe.w,spline2); basis->computeBasis(fe.u,fe.v,fe.w,spline2);
const size_t nbf = basis->order(0)*basis->order(1)*basis->order(2); const size_t nbf = basis->order(0)*basis->order(1)*basis->order(2);
dNdu.resize(nbf,3);
d2Ndu2.resize(nbf,3,3); d2Ndu2.resize(nbf,3,3);
for (size_t n = 1; n <= nbf; n++) { for (size_t n = 1; n <= nbf; n++) {
dNdu(n,1) = spline2.basisDerivs_u[n-1];
dNdu(n,2) = spline2.basisDerivs_v[n-1];
dNdu(n,3) = spline2.basisDerivs_w[n-1];
d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1]; d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1];
d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1]; d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1];
d2Ndu2(n,1,3) = d2Ndu2(n,3,1) = spline2.basisDerivs_uw[n-1]; d2Ndu2(n,1,3) = d2Ndu2(n,3,1) = spline2.basisDerivs_uw[n-1];
@ -298,10 +270,9 @@ bool SplineField3D::hessianFE (const FiniteElement& fe, Matrix& H) const
d2Ndu2(n,3,3) = spline2.basisDerivs_ww[n-1]; d2Ndu2(n,3,3) = spline2.basisDerivs_ww[n-1];
} }
ip.clear();
ASMs3D::scatterInd(basis->numCoefs(0),basis->numCoefs(1),basis->numCoefs(2), ASMs3D::scatterInd(basis->numCoefs(0),basis->numCoefs(1),basis->numCoefs(2),
basis->order(0),basis->order(1),basis->order(2), basis->order(0),basis->order(1),basis->order(2),
spline.left_idx,ip); spline2.left_idx,ip);
} }
Vector Vnod; Vector Vnod;

3
src/ASM/SplineField3D.h
View File

@ -77,6 +77,9 @@ public:
protected: protected:
const Go::SplineVolume* basis; //!< Spline basis description const Go::SplineVolume* basis; //!< Spline basis description
const Go::SplineVolume* vol; //!< Spline geometry description const Go::SplineVolume* vol; //!< Spline geometry description
private:
unsigned char nsd; //!< Number of space dimensions
}; };
#endif #endif

61
src/ASM/SplineFields2D.C
View File

@ -11,6 +11,8 @@
//! //!
//============================================================================== //==============================================================================
#include "GoTools/geometry/SplineSurface.h"
#include "SplineFields2D.h" #include "SplineFields2D.h"
#include "ASMs2D.h" #include "ASMs2D.h"
#include "FiniteElement.h" #include "FiniteElement.h"
@ -18,8 +20,6 @@
#include "Utilities.h" #include "Utilities.h"
#include "Vec3.h" #include "Vec3.h"
#include "GoTools/geometry/SplineSurface.h"
SplineFields2D::SplineFields2D (const ASMs2D* patch, SplineFields2D::SplineFields2D (const ASMs2D* patch,
const RealArray& v, char nbasis, const RealArray& v, char nbasis,
@ -34,6 +34,8 @@ SplineFields2D::SplineFields2D (const ASMs2D* patch,
const int p2 = basis->order_v(); const int p2 = basis->order_v();
nelm = (n1-p1+1)*(n2-p2+1); nelm = (n1-p1+1)*(n2-p2+1);
nsd = patch->getNoSpaceDim();
size_t ofs = 0; size_t ofs = 0;
for (char i = 1; i < nbasis; ++i) for (char i = 1; i < nbasis; ++i)
ofs += patch->getNoNodes(i)*patch->getNoFields(i); ofs += patch->getNoNodes(i)*patch->getNoFields(i);
@ -61,8 +63,9 @@ SplineFields2D::SplineFields2D (const Go::SplineSurface* srf,
const RealArray& v, int cmp, const char* name) const RealArray& v, int cmp, const char* name)
: Fields(name), basis(srf), surf(srf) : Fields(name), basis(srf), surf(srf)
{ {
values = v; nsd = 2; // That is, this constructor can not be used for shells (nsd=3)
nf = cmp; nf = cmp;
values = v;
} }
@ -113,7 +116,7 @@ bool SplineFields2D::valueCoor (const Vec4& x, Vector& vals) const
pt[1] = x[1]; pt[1] = x[1];
pt[2] = x[2]; pt[2] = x[2];
double clo_u, clo_v, dist; double clo_u, clo_v, dist;
#pragma omp critical #pragma omp critical
surf->closestPoint(pt, clo_u, clo_v, clopt, dist, 1e-5); surf->closestPoint(pt, clo_u, clo_v, clopt, dist, 1e-5);
fe.u = clo_u; fe.u = clo_u;
@ -150,9 +153,9 @@ bool SplineFields2D::gradFE (const FiniteElement& fe, Matrix& grad) const
uorder,vorder,spline.left_idx,ip); uorder,vorder,spline.left_idx,ip);
// Evaluate the Jacobian inverse // Evaluate the Jacobian inverse
Matrix Xnod, Jac; Matrix Xnod(nsd,ip.size()), Jac;
Vector Xctrl(&(*surf->coefs_begin()),surf->coefs_end()-surf->coefs_begin()); for (size_t i = 0; i < ip.size(); i++)
utl::gather(ip,surf->dimension(),Xctrl,Xnod); Xnod.fillColumn(1+i,&(*surf->coefs_begin())+surf->dimension()*ip[i]);
utl::Jacobian(Jac,dNdX,Xnod,dNdu); utl::Jacobian(Jac,dNdX,Xnod,dNdu);
// Evaluate the gradient of the solution field at the given point // Evaluate the gradient of the solution field at the given point
@ -189,72 +192,42 @@ bool SplineFields2D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
if (!basis) return false; if (!basis) return false;
if (!surf) return false; if (!surf) return false;
// Order of basis
const int uorder = surf->order_u();
const int vorder = surf->order_v();
const size_t nen = uorder*vorder;
// Evaluate the basis functions at the given point
Go::BasisDerivsSf spline;
Go::BasisDerivsSf2 spline2; Go::BasisDerivsSf2 spline2;
Matrix3D d2Ndu2; Matrix3D d2Ndu2;
Matrix dNdu(nen,2), dNdX;
IntVec ip; IntVec ip;
if (surf == basis) { if (surf == basis) {
#pragma omp critical #pragma omp critical
surf->computeBasis(fe.u,fe.v,spline2); surf->computeBasis(fe.u,fe.v,spline2);
const size_t nen = surf->order_u()*surf->order_v();
d2Ndu2.resize(nen,2,2); d2Ndu2.resize(nen,2,2);
for (size_t n = 1; n <= nen; n++) { for (size_t n = 1; n <= nen; n++) {
dNdu(n,1) = spline2.basisDerivs_u[n-1];
dNdu(n,2) = spline2.basisDerivs_v[n-1];
d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1]; d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1];
d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1]; d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1];
d2Ndu2(n,2,2) = spline2.basisDerivs_vv[n-1]; d2Ndu2(n,2,2) = spline2.basisDerivs_vv[n-1];
} }
ASMs2D::scatterInd(surf->numCoefs_u(),surf->numCoefs_v(), ASMs2D::scatterInd(surf->numCoefs_u(),surf->numCoefs_v(),
uorder,vorder,spline2.left_idx,ip); surf->order_u(),surf->order_v(),
spline2.left_idx,ip);
} }
else { else {
#pragma omp critical
surf->computeBasis(fe.u,fe.v,spline);
for (size_t n = 1; n <= nen; n++) {
dNdu(n,1) = spline.basisDerivs_u[n-1];
dNdu(n,2) = spline.basisDerivs_v[n-1];
}
ASMs2D::scatterInd(surf->numCoefs_u(),surf->numCoefs_v(),
uorder,vorder,spline.left_idx,ip);
}
// Evaluate the Jacobian inverse
Matrix Xnod, Jac;
Vector Xctrl(&(*surf->coefs_begin()),surf->coefs_end()-surf->coefs_begin());
utl::gather(ip,surf->dimension(),Xctrl,Xnod);
utl::Jacobian(Jac,dNdX,Xnod,dNdu);
// Evaluate the gradient of the solution field at the given point
if (basis != surf)
{
// Mixed formulation, the solution uses a different basis than the geometry // Mixed formulation, the solution uses a different basis than the geometry
#pragma omp critical #pragma omp critical
basis->computeBasis(fe.u,fe.v,spline2); basis->computeBasis(fe.u,fe.v,spline2);
const size_t nbf = basis->order_u()*basis->order_v(); const size_t nbf = basis->order_u()*basis->order_v();
dNdu.resize(nbf,2);
d2Ndu2.resize(nbf,2,2); d2Ndu2.resize(nbf,2,2);
for (size_t n = 1; n <= nbf; n++) { for (size_t n = 1; n <= nbf; n++) {
dNdu(n,1) = spline2.basisDerivs_u[n-1];
dNdu(n,2) = spline2.basisDerivs_v[n-1];
d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1]; d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1];
d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1]; d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1];
d2Ndu2(n,2,2) = spline2.basisDerivs_vv[n-1]; d2Ndu2(n,2,2) = spline2.basisDerivs_vv[n-1];
} }
ip.clear();
ASMs2D::scatterInd(basis->numCoefs_u(),basis->numCoefs_v(), ASMs2D::scatterInd(basis->numCoefs_u(),basis->numCoefs_v(),
basis->order_u(),basis->order_v(), basis->order_u(),basis->order_v(),
spline2.left_idx,ip); spline2.left_idx,ip);
} }
Matrix Vnod; Matrix Vnod;

11
src/ASM/SplineFields2D.h
View File

@ -59,22 +59,22 @@ public:
//! \brief Computes the value in a given node/control point. //! \brief Computes the value in a given node/control point.
//! \param[in] node Node number //! \param[in] node Node number
//! \param[out] vals Node values //! \param[out] vals Node values
bool valueNode(size_t node, Vector& vals) const; virtual bool valueNode(size_t node, Vector& vals) const;
//! \brief Computes the value at a given local coordinate. //! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition //! \param[in] fe Finite element definition
//! \param[out] vals Values in local point in given element //! \param[out] vals Values in local point in given element
bool valueFE(const FiniteElement& fe, Vector& vals) const; virtual bool valueFE(const FiniteElement& fe, Vector& vals) const;
//! \brief Computes the value at a given global coordinate. //! \brief Computes the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point //! \param[in] x Global/physical coordinate for point
//! \param[out] vals Values in given physical coordinate //! \param[out] vals Values in given physical coordinate
bool valueCoor(const Vec4& x, Vector& vals) const; virtual bool valueCoor(const Vec4& x, Vector& vals) const;
//! \brief Computes the gradient for a given local coordinate. //! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element //! \param[in] fe Finite element
//! \param[out] grad Gradient of solution in a given local coordinate //! \param[out] grad Gradient of solution in a given local coordinate
bool gradFE(const FiniteElement& fe, Matrix& grad) const; virtual bool gradFE(const FiniteElement& fe, Matrix& grad) const;
//! \brief Computes the hessian for a given local coordinate. //! \brief Computes the hessian for a given local coordinate.
//! \param[in] fe Finite element quantities //! \param[in] fe Finite element quantities
@ -84,6 +84,9 @@ public:
protected: protected:
const Go::SplineSurface* basis; //!< Spline basis description const Go::SplineSurface* basis; //!< Spline basis description
const Go::SplineSurface* surf; //!< Spline geometry description const Go::SplineSurface* surf; //!< Spline geometry description
private:
unsigned char nsd; //!< Number of space dimensions
}; };
#endif #endif

28
src/ASM/SplineFields2Dmx.C
View File

@ -32,8 +32,8 @@ SplineFields2Dmx::SplineFields2Dmx (const ASMs2Dmx* patch,
for (int i = 1; i < *bases.begin(); ++i) for (int i = 1; i < *bases.begin(); ++i)
ofs += patch->getNoNodes(i)*patch->getNoFields(i); ofs += patch->getNoNodes(i)*patch->getNoFields(i);
vit += ofs; vit += ofs;
for (const auto& it : bases) { for (int b : bases) {
size_t nno = patch->getNoNodes(it)*patch->getNoFields(it); size_t nno = patch->getNoNodes(b)*patch->getNoFields(b);
RealArray::const_iterator end = v.size() > nno+ofs ? vit+nno : v.end(); RealArray::const_iterator end = v.size() > nno+ofs ? vit+nno : v.end();
std::copy(vit,end,std::back_inserter(values)); std::copy(vit,end,std::back_inserter(values));
vit += nno; vit += nno;
@ -58,8 +58,8 @@ bool SplineFields2Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
// Evaluate the basis functions at the given point // Evaluate the basis functions at the given point
auto vit = values.begin(); auto vit = values.begin();
auto rit = vals.begin(); auto rit = vals.begin();
for (const auto& it : bases) { for (int b : bases) {
Go::SplineSurface* basis = surf->getBasis(it); Go::SplineSurface* basis = surf->getBasis(b);
Go::BasisPtsSf spline; Go::BasisPtsSf spline;
#pragma omp critical #pragma omp critical
basis->computeBasis(fe.u,fe.v,spline); basis->computeBasis(fe.u,fe.v,spline);
@ -71,11 +71,11 @@ bool SplineFields2Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
spline.left_idx,ip); spline.left_idx,ip);
Matrix Vnod; Matrix Vnod;
utl::gather(ip,1,Vector(&*vit,surf->getNoNodes(it)),Vnod); utl::gather(ip,1,Vector(&*vit,surf->getNoNodes(b)),Vnod);
Vector val2; Vector val2;
Vnod.multiply(spline.basisValues,val2); // vals = Vnod * basisValues Vnod.multiply(spline.basisValues,val2); // vals = Vnod * basisValues
*rit++ = val2.front(); *rit++ = val2.front();
vit += surf->getNoNodes(it); vit += surf->getNoNodes(b);
} }
return true; return true;
@ -108,16 +108,16 @@ bool SplineFields2Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
uorder,vorder,spline.left_idx,ip); uorder,vorder,spline.left_idx,ip);
// Evaluate the Jacobian inverse // Evaluate the Jacobian inverse
Matrix Xnod, Jac; Matrix Xnod(surf->getNoSpaceDim(),ip.size()), Jac;
Vector Xctrl(&(*gsurf->coefs_begin()),gsurf->coefs_end()-gsurf->coefs_begin()); for (size_t i = 0; i < ip.size(); i++)
utl::gather(ip,gsurf->dimension(),Xctrl,Xnod); Xnod.fillColumn(1+i,&(*gsurf->coefs_begin())+gsurf->dimension()*ip[i]);
utl::Jacobian(Jac,dNdX,Xnod,dNdu); utl::Jacobian(Jac,dNdX,Xnod,dNdu);
// Evaluate the gradient of the solution field at the given point // Evaluate the gradient of the solution field at the given point
auto vit = values.begin(); auto vit = values.begin();
size_t row = 1; size_t row = 1;
for (const auto& it : bases) { for (int b : bases) {
const Go::SplineSurface* basis = surf->getBasis(it); const Go::SplineSurface* basis = surf->getBasis(b);
#pragma omp critical #pragma omp critical
basis->computeBasis(fe.u,fe.v,spline); basis->computeBasis(fe.u,fe.v,spline);
@ -135,13 +135,13 @@ bool SplineFields2Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
basis->order_u(),basis->order_v(), basis->order_u(),basis->order_v(),
spline.left_idx,ip); spline.left_idx,ip);
utl::gather(ip,1,Vector(&*vit, surf->getNoNodes(it)* const size_t nval = surf->getNoNodes(b)*surf->getNoFields(b);
surf->getNoFields(it)),Xnod); utl::gather(ip,1,Vector(&*vit,nval),Xnod);
Matrix grad2; Matrix grad2;
grad2.multiply(Xnod,dNdX); // grad = Xnod * dNdX grad2.multiply(Xnod,dNdX); // grad = Xnod * dNdX
grad(row,1) = grad2(1,1); grad(row,1) = grad2(1,1);
grad(row++,2) = grad2(1,2); grad(row++,2) = grad2(1,2);
vit += surf->getNoNodes(it)*surf->getNoFields(it); vit += nval;
} }
return true; return true;

57
src/ASM/SplineFields3D.C
View File

@ -36,6 +36,8 @@ SplineFields3D::SplineFields3D (const ASMs3D* patch,
const int p3 = basis->order(2); const int p3 = basis->order(2);
nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1); nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1);
nsd = patch->getNoSpaceDim();
size_t ofs = 0; size_t ofs = 0;
for (char i = 1; i < nbasis; ++i) for (char i = 1; i < nbasis; ++i)
ofs += patch->getNoNodes(i)*patch->getNoFields(i); ofs += patch->getNoNodes(i)*patch->getNoFields(i);
@ -62,8 +64,9 @@ SplineFields3D::SplineFields3D (const Go::SplineVolume* svol,
const RealArray& v, int cmp, const char* name) const RealArray& v, int cmp, const char* name)
: Fields(name), basis(svol), vol(svol) : Fields(name), basis(svol), vol(svol)
{ {
values = v; nsd = 3;
nf = cmp; nf = cmp;
values = v;
} }
@ -115,7 +118,7 @@ bool SplineFields3D::valueCoor (const Vec4& x, Vector& vals) const
pt[1] = x[1]; pt[1] = x[1];
pt[2] = x[2]; pt[2] = x[2];
double clo_u, clo_v, clo_w, dist; double clo_u, clo_v, clo_w, dist;
#pragma omp critical #pragma omp critical
vol->closestPoint(pt, clo_u, clo_v, clo_w, clopt, dist, 1e-5); vol->closestPoint(pt, clo_u, clo_v, clo_w, clopt, dist, 1e-5);
fe.u = clo_u; fe.u = clo_u;
@ -155,9 +158,9 @@ bool SplineFields3D::gradFE (const FiniteElement& fe, Matrix& grad) const
uorder,vorder,worder,spline.left_idx,ip); uorder,vorder,worder,spline.left_idx,ip);
// Evaluate the Jacobian inverse // Evaluate the Jacobian inverse
Matrix Xnod, Jac; Matrix Xnod(nsd,ip.size()), Jac;
Vector Xctrl(&(*vol->coefs_begin()),vol->coefs_end()-vol->coefs_begin()); for (size_t i = 0; i < ip.size(); i++)
utl::gather(ip,vol->dimension(),Xctrl,Xnod); Xnod.fillColumn(1+i,&(*vol->coefs_begin())+vol->dimension()*ip[i]);
utl::Jacobian(Jac,dNdX,Xnod,dNdu); utl::Jacobian(Jac,dNdX,Xnod,dNdu);
// Evaluate the gradient of the solution field at the given point // Evaluate the gradient of the solution field at the given point
@ -195,27 +198,17 @@ bool SplineFields3D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
if (!basis) return false; if (!basis) return false;
if (!vol) return false; if (!vol) return false;
const int uorder = vol->order(0);
const int vorder = vol->order(1);
const int worder = vol->order(2);
const size_t nen = uorder*vorder*worder;
// Evaluate the basis functions at the given point // Evaluate the basis functions at the given point
Go::BasisDerivs spline;
Go::BasisDerivs2 spline2; Go::BasisDerivs2 spline2;
Matrix3D d2Ndu2; Matrix3D d2Ndu2;
Matrix dNdu, dNdX;
IntVec ip; IntVec ip;
if (vol == basis) { if (vol == basis) {
#pragma omp critical #pragma omp critical
vol->computeBasis(fe.u,fe.v,fe.w,spline2); vol->computeBasis(fe.u,fe.v,fe.w,spline2);
dNdu.resize(nen,3); const size_t nen = vol->order(0)*vol->order(1)*vol->order(2);
d2Ndu2.resize(nen,3,3); d2Ndu2.resize(nen,3,3);
for (size_t n = 1; n <= nen; n++) { for (size_t n = 1; n <= nen; n++) {
dNdu(n,1) = spline2.basisDerivs_u[n-1];
dNdu(n,2) = spline2.basisDerivs_v[n-1];
dNdu(n,3) = spline2.basisDerivs_w[n-1];
d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1]; d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1];
d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1]; d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1];
d2Ndu2(n,1,3) = d2Ndu2(n,3,1) = spline2.basisDerivs_uw[n-1]; d2Ndu2(n,1,3) = d2Ndu2(n,3,1) = spline2.basisDerivs_uw[n-1];
@ -225,42 +218,17 @@ bool SplineFields3D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
} }
ASMs3D::scatterInd(vol->numCoefs(0),vol->numCoefs(1),vol->numCoefs(2), ASMs3D::scatterInd(vol->numCoefs(0),vol->numCoefs(1),vol->numCoefs(2),
uorder,vorder,worder,spline2.left_idx,ip); vol->order(0),vol->order(1),vol->order(2),
spline2.left_idx,ip);
} }
else { else {
#pragma omp critical
vol->computeBasis(fe.u,fe.v,fe.w,spline);
dNdu.resize(nen,3);
for (size_t n = 1; n <= nen; n++) {
dNdu(n,1) = spline.basisDerivs_u[n-1];
dNdu(n,2) = spline.basisDerivs_v[n-1];
dNdu(n,3) = spline.basisDerivs_w[n-1];
}
ASMs3D::scatterInd(vol->numCoefs(0),vol->numCoefs(1),vol->numCoefs(2),
uorder,vorder,worder,spline.left_idx,ip);
}
// Evaluate the Jacobian inverse
Matrix Xnod, Jac;
Vector Xctrl(&(*vol->coefs_begin()),vol->coefs_end()-vol->coefs_begin());
utl::gather(ip,vol->dimension(),Xctrl,Xnod);
utl::Jacobian(Jac,dNdX,Xnod,dNdu);
// Evaluate the gradient of the solution field at the given point
if (basis != vol) {
// Mixed formulation, the solution uses a different basis than the geometry // Mixed formulation, the solution uses a different basis than the geometry
#pragma omp critical #pragma omp critical
basis->computeBasis(fe.u,fe.v,fe.w,spline2); basis->computeBasis(fe.u,fe.v,fe.w,spline2);
const size_t nbf = basis->order(0)*basis->order(1)*basis->order(2); const size_t nbf = basis->order(0)*basis->order(1)*basis->order(2);
dNdu.resize(nbf,3);
d2Ndu2.resize(nbf,3,3); d2Ndu2.resize(nbf,3,3);
for (size_t n = 1; n <= nbf; n++) { for (size_t n = 1; n <= nbf; n++) {
dNdu(n,1) = spline2.basisDerivs_u[n-1];
dNdu(n,2) = spline2.basisDerivs_v[n-1];
dNdu(n,3) = spline2.basisDerivs_w[n-1];
d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1]; d2Ndu2(n,1,1) = spline2.basisDerivs_uu[n-1];
d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1]; d2Ndu2(n,1,2) = d2Ndu2(n,2,1) = spline2.basisDerivs_uv[n-1];
d2Ndu2(n,1,3) = d2Ndu2(n,3,1) = spline2.basisDerivs_uw[n-1]; d2Ndu2(n,1,3) = d2Ndu2(n,3,1) = spline2.basisDerivs_uw[n-1];
@ -269,10 +237,9 @@ bool SplineFields3D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
d2Ndu2(n,3,3) = spline2.basisDerivs_ww[n-1]; d2Ndu2(n,3,3) = spline2.basisDerivs_ww[n-1];
} }
ip.clear();
ASMs3D::scatterInd(basis->numCoefs(0),basis->numCoefs(1),basis->numCoefs(2), ASMs3D::scatterInd(basis->numCoefs(0),basis->numCoefs(1),basis->numCoefs(2),
basis->order(0),basis->order(1),basis->order(2), basis->order(0),basis->order(1),basis->order(2),
spline.left_idx,ip); spline2.left_idx,ip);
} }
Matrix Vnod; Matrix Vnod;

11
src/ASM/SplineFields3D.h
View File

@ -59,22 +59,22 @@ public:
//! \brief Computes the value in a given node/control point. //! \brief Computes the value in a given node/control point.
//! \param[in] node Node number //! \param[in] node Node number
//! \param[out] vals Node values //! \param[out] vals Node values
bool valueNode(size_t node, Vector& vals) const; virtual bool valueNode(size_t node, Vector& vals) const;
//! \brief Computes the value at a given local coordinate. //! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition //! \param[in] fe Finite element definition
//! \param[out] vals Values in local point in given element //! \param[out] vals Values in local point in given element
bool valueFE(const FiniteElement& fe, Vector& vals) const; virtual bool valueFE(const FiniteElement& fe, Vector& vals) const;
//! \brief Computes the value at a given global coordinate. //! \brief Computes the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point //! \param[in] x Global/physical coordinate for point
//! \param[out] vals Values in given physical coordinate //! \param[out] vals Values in given physical coordinate
bool valueCoor(const Vec4& x, Vector& vals) const; virtual bool valueCoor(const Vec4& x, Vector& vals) const;
//! \brief Computes the gradient for a given local coordinate. //! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element //! \param[in] fe Finite element
//! \param[out] grad Gradient of solution in a given local coordinate //! \param[out] grad Gradient of solution in a given local coordinate
bool gradFE(const FiniteElement& fe, Matrix& grad) const; virtual bool gradFE(const FiniteElement& fe, Matrix& grad) const;
//! \brief Computes the hessian for a given local coordinate. //! \brief Computes the hessian for a given local coordinate.
//! \param[in] fe Finite element quantities //! \param[in] fe Finite element quantities
@ -84,6 +84,9 @@ public:
protected: protected:
const Go::SplineVolume* basis; //!< Spline basis description const Go::SplineVolume* basis; //!< Spline basis description
const Go::SplineVolume* vol; //!< Spline geometry description const Go::SplineVolume* vol; //!< Spline geometry description
private:
unsigned char nsd; //!< Number of space dimensions
}; };
#endif #endif

30
src/ASM/SplineFields3Dmx.C
View File

@ -32,8 +32,8 @@ SplineFields3Dmx::SplineFields3Dmx (const ASMs3Dmx* patch,
for (int i = 1; i < *bases.begin(); ++i) for (int i = 1; i < *bases.begin(); ++i)
ofs += patch->getNoNodes(i)*patch->getNoFields(i); ofs += patch->getNoNodes(i)*patch->getNoFields(i);
vit += ofs; vit += ofs;
for (const auto& it : bases) { for (int b : bases) {
size_t nno = patch->getNoFields(it)*patch->getNoNodes(it); size_t nno = patch->getNoFields(b)*patch->getNoNodes(b);
RealArray::const_iterator end = v.size() > nno+ofs ? vit+nno : v.end(); RealArray::const_iterator end = v.size() > nno+ofs ? vit+nno : v.end();
std::copy(vit,end,std::back_inserter(values)); std::copy(vit,end,std::back_inserter(values));
vit += nno; vit += nno;
@ -58,8 +58,8 @@ bool SplineFields3Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
// Evaluate the basis functions at the given point // Evaluate the basis functions at the given point
auto vit = values.begin(); auto vit = values.begin();
auto rit = vals.begin(); auto rit = vals.begin();
for (const auto& it : bases) { for (int b : bases) {
Go::SplineVolume* basis = svol->getBasis(it); Go::SplineVolume* basis = svol->getBasis(b);
Go::BasisPts spline; Go::BasisPts spline;
#pragma omp critical #pragma omp critical
basis->computeBasis(fe.u,fe.v,fe.w,spline); basis->computeBasis(fe.u,fe.v,fe.w,spline);
@ -72,11 +72,11 @@ bool SplineFields3Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
spline.left_idx,ip); spline.left_idx,ip);
Matrix Vnod; Matrix Vnod;
utl::gather(ip,1,Vector(&*vit,svol->getNoNodes(it)),Vnod); utl::gather(ip,1,Vector(&*vit,svol->getNoNodes(b)),Vnod);
Vector val2; Vector val2;
Vnod.multiply(spline.basisValues,val2); // vals = Vnod * basisValues Vnod.multiply(spline.basisValues,val2); // vals = Vnod * basisValues
*rit++ = val2.front(); *rit++ = val2.front();
vit += svol->getNoNodes(it); vit += svol->getNoNodes(b);
} }
return true; return true;
@ -111,16 +111,16 @@ bool SplineFields3Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
uorder,vorder,worder,spline.left_idx,ip); uorder,vorder,worder,spline.left_idx,ip);
// Evaluate the Jacobian inverse // Evaluate the Jacobian inverse
Matrix Xnod, Jac; Matrix Xnod(svol->getNoSpaceDim(),ip.size()), Jac;
Vector Xctrl(&(*gvol->coefs_begin()),gvol->coefs_end()-gvol->coefs_begin()); for (size_t i = 0; i < ip.size(); i++)
utl::gather(ip,gvol->dimension(),Xctrl,Xnod); Xnod.fillColumn(1+i,&(*gvol->coefs_begin())+gvol->dimension()*ip[i]);
utl::Jacobian(Jac,dNdX,Xnod,dNdu); utl::Jacobian(Jac,dNdX,Xnod,dNdu);
// Evaluate the gradient of the solution field at the given point // Evaluate the gradient of the solution field at the given point
auto vit = values.begin(); auto vit = values.begin();
size_t row = 1; size_t row = 1;
for (const auto& it : bases) { for (int b : bases) {
const Go::SplineVolume* basis = svol->getBasis(it); const Go::SplineVolume* basis = svol->getBasis(b);
#pragma omp critical #pragma omp critical
basis->computeBasis(fe.u,fe.v,fe.w,spline); basis->computeBasis(fe.u,fe.v,fe.w,spline);
@ -139,14 +139,14 @@ bool SplineFields3Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
basis->order(0),basis->order(1),basis->order(2), basis->order(0),basis->order(1),basis->order(2),
spline.left_idx,ip); spline.left_idx,ip);
utl::gather(ip,1,Vector(&*vit, svol->getNoNodes(it)* const size_t nval = svol->getNoNodes(b)*svol->getNoFields(b);
svol->getNoFields(it)),Xnod); utl::gather(ip,1,Vector(&*vit,nval),Xnod);
Matrix grad2; Matrix grad2;
grad2.multiply(Xnod,dNdX); // grad = Xnod * dNdX grad2.multiply(Xnod,dNdX); // grad = Xnod * dNdX
grad(row,1) = grad2(1,1); grad(row,1) = grad2(1,1);
grad(row,2) = grad(1,2); grad(row,2) = grad(1,2);
grad(row++,3) = grad(1,3); grad(row++,3) = grad(1,3);
vit += svol->getNoNodes(it)*svol->getNoFields(it); vit += nval;
} }
return true; return true;