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

changed: use BasisFunctionValues also in boundary integrals

Browse Source 
this allows simplifying the MxFiniteElement::(Hessian|Jacobian) methods.

add a type alias and a helper class to make it convenient.
This commit is contained in:
Arne Morten Kvarving 2023-09-15 15:02:32 +02:00
parent 1e96b51374
commit f131400ec8
10 changed files with 191 additions and 150 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/ASM/ASMs2D.C
View File

@ -1423,7 +1423,8 @@ bool ASMs2D::getElementCoordinatesPrm (Matrix& X, double u, double v) const
}
#if SP_DEBUG > 2
std::cout <<"\nCoordinates for element "<< iel << X << std::endl;
std::cout <<"\nCoordinates for element containing parameters ("
<< u <<"," << v << "):" << X << std::endl;
#endif
return true;
}

40
src/ASM/ASMs2Dmx.C
View File

@ -551,7 +551,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
fe.v = param[1] = cache.getParam(1,i2-p2,j);
// Fetch basis function derivatives at current integration point
std::vector<const BasisFunctionVals*> bfs(myCache.size());
BasisValuesPtrs bfs(myCache.size());
for (size_t b = 0; b < myCache.size(); ++b) {
bfs[b] = &myCache[b]->getVals(iel-1, ip);
if (b < m_basis.size())
@ -560,11 +560,11 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,Xnod,itgBasis,&bfs))
if (!fe.Jacobian(Jac,Xnod,itgBasis,bfs))
continue; // skip singular points
// Compute Hessian of coordinate mapping and 2nd order derivatives
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,itgBasis,&bfs))
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,itgBasis,bfs))
ok = false;
// Compute G-matrix
@ -662,8 +662,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
fe.v = gpar[1](1,1);
double param[3] = { fe.u, fe.v, 0.0 };
Matrices dNxdu(splinex.size());
Vector Ng;
BasisValues bfs(splinex.size());
Matrix Xnod, Jac;
Vec4 X(param,time.t);
Vec3 normal;
@ -729,21 +728,21 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
}
if (separateGeometry)
SplineUtils::extractBasis(splinex.back()[ip],Ng,dNxdu.back());
SplineUtils::extractBasis(splinex.back()[ip],bfs.back().N,bfs.back().dNdu);
// Fetch basis function derivatives at current integration point
for (size_t b = 0; b < m_basis.size(); ++b)
SplineUtils::extractBasis(splinex[b][ip],fe.basis(b+1),dNxdu[b]);
SplineUtils::extractBasis(splinex[b][ip],fe.basis(b+1),bfs[b].dNdu);
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,dNxdu,t1,t2))
if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,bfs,t1,t2))
continue; // skip singular points
if (edgeDir < 0) normal *= -1.0;
// Cartesian coordinates of current integration point
X.assign(Xnod * (separateGeometry ? Ng : fe.basis(itgBasis)));
X.assign(Xnod * (separateGeometry ? bfs.back().N : fe.basis(itgBasis)));
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= dS*wg[i];
@ -811,8 +810,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
std::copy(elem_size.begin(),elem_size.end(),std::back_inserter(elem_sizes2));
MxFiniteElement fe(elem_sizes2);
Matrix dNdu, Xnod, Jac;
Vector dN;
Matrix Xnod, Jac;
Vec4 X(nullptr,time.t);
Vec3 normal;
double u[2], v[2];
@ -906,17 +904,17 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
}
// Fetch basis function derivatives at current integration point
Matrices dNxdu(2*nB);
BasisValues bfs(2*nB);
for (size_t b = 1; b <= nB; b++) {
Go::BasisDerivsSf spline;
this->getBasis(b)->computeBasis(fe.u, fe.v, spline, edgeDir < 0);
SplineUtils::extractBasis(spline, fe.basis(b), dNxdu[b-1]);
this->getBasis(b)->computeBasis(fe.u, fe.v, spline, edgeDir > 0);
SplineUtils::extractBasis(spline, fe.basis(nB+b), dNxdu[nB+b-1]);
this->getBasis(b)->computeBasis(fe.u,fe.v,spline,edgeDir < 0);
SplineUtils::extractBasis(spline,fe.basis(b),bfs[b-1].dNdu);
this->getBasis(b)->computeBasis(fe.u,fe.v,spline,edgeDir > 0);
SplineUtils::extractBasis(spline,fe.basis(nB+b),bfs[nB+b-1].dNdu);
}
// Compute basis function derivatives and the edge normal
if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,dNxdu,t1,t2,nB))
if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,bfs,t1,t2,nB))
continue; // skip singular points
if (edgeDir < 0) normal *= -1.0;
@ -1056,7 +1054,7 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
MxFiniteElement fe(elem_size,firstIp);
Vector solPt;
Matrices dNxdu(m_basis.size());
BasisValues bfs(m_basis.size());
Matrix Jac;
// Evaluate the secondary solution field at each point
@ -1073,7 +1071,7 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
splinex[b][i].left_idx,ip[b]);
// Fetch associated control point coordinates
if (b == (size_t)itgBasis-1)
if (b == static_cast<size_t>(itgBasis)-1)
utl::gather(ip[itgBasis-1], nsd, Xnod, Xtmp);
for (int& c : ip[b]) c += ofs;
@ -1086,11 +1084,11 @@ bool ASMs2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
// Fetch basis function derivatives at current integration point
for (size_t b = 0; b < m_basis.size(); ++b)
SplineUtils::extractBasis(splinex[b][i],fe.basis(b+1),dNxdu[b]);
SplineUtils::extractBasis(splinex[b][i],fe.basis(b+1),bfs[b].dNdu);
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,Xtmp,itgBasis,nullptr,&dNxdu))
if (!fe.Jacobian(Jac,Xtmp,itgBasis,bfs))
continue; // skip singular points
// Cartesian coordinates of current integration point

16
src/ASM/ASMs2DmxLag.C
View File

@ -312,14 +312,14 @@ bool ASMs2DmxLag::integrate (Integrand& integrand,
fe.eta = xg[1][j];
// Fetch basis function derivatives at current integration point
std::vector<const BasisFunctionVals*> bfs(this->getNoBasis());
BasisValuesPtrs bfs(this->getNoBasis());
for (size_t b = 0; b < this->getNoBasis(); ++b) {
bfs[b] = &myCache[b]->getVals(iel-1,ip);
fe.basis(b+1) = bfs[b]->N;
}
// Compute Jacobian inverse of coordinate mapping and derivatives
if (!fe.Jacobian(Jac,Xnod,itgBasis,&bfs))
if (!fe.Jacobian(Jac,Xnod,itgBasis,bfs))
continue; // skip singular points
// Cartesian coordinates of current integration point
@ -377,7 +377,7 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, int lIndex,
size_t firstp = iit == firstBp.end() ? 0 : iit->second;
MxFiniteElement fe(elem_size);
Matrices dNxdu(nxx.size());
BasisValues bfs(nxx.size());
Matrix Xnod, Jac;
Vec4 X(nullptr,time.t);
Vec3 normal;
@ -422,13 +422,13 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, int lIndex,
// Compute the basis functions and their derivatives, using
// tensor product of one-dimensional Lagrange polynomials
for (size_t b = 0; b < nxx.size(); ++b)
if (!Lagrange::computeBasis(fe.basis(b+1),dNxdu[b],
if (!Lagrange::computeBasis(fe.basis(b+1),bfs[b].dNdu,
elem_sizes[b][0],xi[0],
elem_sizes[b][1],xi[1]))
ok = false;
// Compute basis function derivatives and the edge normal
if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,dNxdu,t1,t2))
if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,bfs,t1,t2))
continue; // skip singular points
if (edgeDir < 0) normal *= -1.0;
@ -499,7 +499,7 @@ bool ASMs2DmxLag::evalSolution (Matrix& sField, const IntegrandBase& integrand,
MxFiniteElement fe(elem_size);
Vector solPt;
Vectors globSolPt(nPoints);
Matrices dNxdu(nxx.size());
BasisValues bfs(nxx.size());
Matrix Xnod, Jac;
// Evaluate the secondary solution field at each point
@ -521,14 +521,14 @@ bool ASMs2DmxLag::evalSolution (Matrix& sField, const IntegrandBase& integrand,
double xi = -1.0 + i*incx;
double eta = -1.0 + j*incy;
for (size_t b = 0; b < nxx.size(); ++b)
if (!Lagrange::computeBasis(fe.basis(b+1),dNxdu[b],
if (!Lagrange::computeBasis(fe.basis(b+1),bfs[b].dNdu,
elem_sizes[b][0],xi,
elem_sizes[b][1],eta))
return false;
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,Xnod,itgBasis,nullptr,&dNxdu))
if (!fe.Jacobian(Jac,Xnod,itgBasis,bfs))
continue; // skip singular points
// Now evaluate the solution field

38
src/ASM/ASMs3Dmx.C
View File

@ -582,7 +582,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
fe.w = param[2] = cache.getParam(2,i3-p3,k);
// Fetch basis function derivatives at current integration point
std::vector<const BasisFunctionVals*> bfs(myCache.size());
BasisValuesPtrs bfs(myCache.size());
for (size_t b = 0; b < myCache.size(); ++b) {
bfs[b] = &myCache[b]->getVals(iel-1, ip);
if (b < m_basis.size())
@ -591,11 +591,11 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,Xnod,itgBasis,&bfs))
if (!fe.Jacobian(Jac,Xnod,itgBasis,bfs))
continue; // skip singular points
// Compute Hessian of coordinate mapping and 2nd order derivatives
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,itgBasis,&bfs))
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,itgBasis,bfs))
ok = false;
// Compute G-matrix
@ -716,8 +716,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
fe.v = gpar[1](1,1);
fe.w = gpar[2](1,1);
Matrices dNxdu(splinex.size());
Vector Ng;
BasisValues bfs(splinex.size());
Matrix Xnod, Jac;
double param[3] = { fe.u, fe.v, fe.w };
Vec4 X(param,time.t);
@ -806,21 +805,21 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
}
if (separateGeometry)
SplineUtils::extractBasis(splinex.back()[ip],Ng,dNxdu.back());
SplineUtils::extractBasis(splinex.back()[ip],bfs.back().N,bfs.back().dNdu);
// Fetch basis function derivatives at current integration point
for (size_t b = 0; b < m_basis.size(); ++b)
SplineUtils::extractBasis(splinex[b][ip],fe.basis(b+1),dNxdu[b]);
SplineUtils::extractBasis(splinex[b][ip],fe.basis(b+1),bfs[b].dNdu);
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,dNxdu,t1,t2))
if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,bfs,t1,t2))
continue; // skip singular points
if (faceDir < 0) normal *= -1.0;
// Cartesian coordinates of current integration point
X.assign(Xnod * (separateGeometry ? Ng : fe.basis(itgBasis)));
X.assign(Xnod * (separateGeometry ? bfs.back().N : fe.basis(itgBasis)));
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= dA*wg[i]*wg[j];
@ -892,8 +891,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
std::copy(elem_size.begin(), elem_size.end(), std::back_inserter(elem_sizes2));
MxFiniteElement fe(elem_sizes2);
Matrix dNdu, Xnod, Jac;
Vector dN;
Matrix Xnod, Jac;
Vec4 X(nullptr,time.t);
Vec3 normal;
double u[2], v[2], w[2];
@ -1021,17 +1019,17 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
}
// Fetch basis function derivatives at current integration point
Matrices dNxdu(2*nB);
BasisValues bfs(2*nB);
for (size_t b = 1; b <= nB; b++) {
Go::BasisDerivs spline;
this->getBasis(b)->computeBasis(fe.u, fe.v, fe.w, spline, faceDir < 0);
SplineUtils::extractBasis(spline, fe.basis(b), dNxdu[b-1]);
this->getBasis(b)->computeBasis(fe.u, fe.v, fe.w, spline, faceDir > 0);
SplineUtils::extractBasis(spline, fe.basis(nB+b), dNxdu[nB+b-1]);
this->getBasis(b)->computeBasis(fe.u,fe.v,fe.w,spline,faceDir < 0);
SplineUtils::extractBasis(spline,fe.basis(b),bfs[b-1].dNdu);
this->getBasis(b)->computeBasis(fe.u,fe.v,fe.w,spline,faceDir > 0);
SplineUtils::extractBasis(spline,fe.basis(nB+b),bfs[nB+b-1].dNdu);
}
// Compute basis function derivatives and the edge normal
if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,dNxdu,t1,t2,nB))
if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,bfs,t1,t2,nB))
continue; // skip singular points
if (faceDir < 0) normal *= -1.0;
@ -1173,7 +1171,7 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
MxFiniteElement fe(elem_size,firstIp);
Vector solPt;
Matrices dNxdu(m_basis.size());
BasisValues bfs(m_basis.size());
Matrix Jac;
// Evaluate the secondary solution field at each point
@ -1204,11 +1202,11 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
// Fetch basis function derivatives at current integration point
for (size_t b = 0; b < m_basis.size(); ++b)
SplineUtils::extractBasis(splinex[b][i],fe.basis(b+1),dNxdu[b]);
SplineUtils::extractBasis(splinex[b][i],fe.basis(b+1),bfs[b].dNdu);
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinate
if (!fe.Jacobian(Jac,Xtmp,itgBasis,nullptr,&dNxdu))
if (!fe.Jacobian(Jac,Xtmp,itgBasis,bfs))
continue; // skip singular points
// Cartesian coordinates of current integration point

17
src/ASM/ASMs3DmxLag.C
View File

@ -18,7 +18,6 @@
#include "GlobalIntegral.h"
#include "LocalIntegral.h"
#include "IntegrandBase.h"
#include "CoordinateMapping.h"
#include "GaussQuadrature.h"
#include <numeric>
@ -348,14 +347,14 @@ bool ASMs3DmxLag::integrate (Integrand& integrand,
fe.zeta = xg[2][k];
// Compute basis function derivatives at current integration point
std::vector<const BasisFunctionVals*> bfs(this->getNoBasis());
BasisValuesPtrs bfs(this->getNoBasis());
for (size_t b = 0; b < this->getNoBasis(); ++b) {
bfs[b] = &myCache[b]->getVals(iel-1,ip);
fe.basis(b+1) = bfs[b]->N;
}
// Compute Jacobian inverse of coordinate mapping and derivatives
if (!fe.Jacobian(Jac,Xnod,itgBasis,&bfs))
if (!fe.Jacobian(Jac,Xnod,itgBasis,bfs))
continue; // skip singular points
// Cartesian coordinates of current integration point
@ -430,7 +429,7 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, int lIndex,
for (size_t t = 0; t < threadGrp[g].size(); t++)
{
MxFiniteElement fe(elem_size);
Matrices dNxdu;
BasisValues bfs(elem_size.size());
Matrix Xnod, Jac;
Vec4 X(nullptr,time.t);
Vec3 normal;
@ -487,14 +486,14 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, int lIndex,
// Compute the basis functions and their derivatives, using
// tensor product of one-dimensional Lagrange polynomials
for (size_t b = 0; b < nxx.size(); ++b)
if (!Lagrange::computeBasis(fe.basis(b+1),dNxdu[b],
if (!Lagrange::computeBasis(fe.basis(b+1),bfs[b].dNdu,
elem_sizes[b][0],xi[0],
elem_sizes[b][1],xi[1],
elem_sizes[b][2],xi[2]))
ok = false;
// Compute basis function derivatives and the edge normal
if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,dNxdu,t1,t2))
if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,bfs,t1,t2))
continue; // skip singular points
if (faceDir < 0) normal *= -1.0;
@ -566,7 +565,7 @@ bool ASMs3DmxLag::evalSolution (Matrix& sField, const IntegrandBase& integrand,
MxFiniteElement fe(elem_size);
Vector solPt;
Vectors globSolPt(nPoints);
Matrices dNxdu(nxx.size());
BasisValues bfs(nxx.size());
Matrix Xnod, Jac;
// Evaluate the secondary solution field at each point
@ -590,7 +589,7 @@ bool ASMs3DmxLag::evalSolution (Matrix& sField, const IntegrandBase& integrand,
fe.eta = -1.0 + j*incy;
fe.zeta = -1.0 + k*incz;
for (size_t b = 0; b < nxx.size(); ++b)
if (!Lagrange::computeBasis(fe.basis(b+1),dNxdu[b],
if (!Lagrange::computeBasis(fe.basis(b+1),bfs[b].dNdu,
elem_sizes[b][0],fe.xi,
elem_sizes[b][1],fe.eta,
elem_sizes[b][2],fe.zeta))
@ -598,7 +597,7 @@ bool ASMs3DmxLag::evalSolution (Matrix& sField, const IntegrandBase& integrand,
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,Xnod,itgBasis,nullptr,&dNxdu))
if (!fe.Jacobian(Jac,Xnod,itgBasis,bfs))
continue; // skip singular points
// Now evaluate the solution field

59
src/ASM/BasisFunctionVals.h
View File

@ -29,4 +29,63 @@ struct BasisFunctionVals
Matrix4D d3Ndu3; //!< Third order derivatives of basis functions
};
//! \brief Convenience type alias
using BasisValuesPtrs = std::vector<const BasisFunctionVals*>;
//! \brief Utility class holding a vector of basis function values.
//! \details This is a helper class that allows us to avoid pointer
//! semantics in code, while allowing to pass as
//! pointers into mapping functions. This is necessary
//! to avoid additional copying in code using the
//! basis function cache.
//! While it overrides most common vector operations,
//! be careful with advanced insertions.
class BasisValues : public std::vector<BasisFunctionVals>
{
public:
//! \brief Constructor resizing to a given size.
BasisValues(size_t size) : std::vector<BasisFunctionVals>(size)
{
pointers.reserve(this->size());
for (const BasisFunctionVals& val : *this)
pointers.push_back(&val);
}
//! \brief Destructor.
//! \details Need to manually implement to call parent dtor.
~BasisValues()
{
this->clear();
}
//! \brief Clears the container.
void clear() noexcept
{
this->std::vector<BasisFunctionVals>::clear();
pointers.clear();
}
//! \brief Push back for reference.
void push_back(const BasisFunctionVals& v)
{
this->std::vector<BasisFunctionVals>::push_back(v);
pointers.push_back(&this->back());
}
//! \brief Push back for r-reference.
void push_back(BasisFunctionVals&& v)
{
this->std::vector<BasisFunctionVals>::push_back(v);
pointers.push_back(&this->back());
}
//! \brief Cast to a vector with pointers to our values.
operator const BasisValuesPtrs& () const { return pointers; }
private:
BasisValuesPtrs pointers; //!< Vector of pointers to our values
};
#endif

36
src/ASM/FiniteElement.C
View File

@ -77,22 +77,20 @@ std::ostream& MxFiniteElement::write (std::ostream& os) const
bool MxFiniteElement::Jacobian (Matrix& Jac, const Matrix& Xnod,
unsigned short int gBasis,
const std::vector<const BasisFunctionVals*>* bf,
const std::vector<Matrix>* dNxdu)
const BasisValuesPtrs& bfs)
{
size_t nBasis = bf ? bf->size() : dNxdu->size();
const bool separateGeometry = nBasis > this->getNoBasis();
const bool separateGeometry = bfs.size() > this->getNoBasis();
if (separateGeometry)
gBasis = nBasis;
gBasis = bfs.size();
detJxW = utl::Jacobian(Jac, this->grad(gBasis), Xnod,
bf ? (*bf)[gBasis-1]->dNdu : (*dNxdu)[gBasis-1],
bfs[gBasis-1]->dNdu,
!separateGeometry);
if (detJxW == 0.0) return false; // singular point
for (size_t b = 1; b <= this->getNoBasis(); b++)
if (b != gBasis || separateGeometry)
this->grad(b).multiply(bf ? (*bf)[b-1]->dNdu : (*dNxdu)[b-1], Jac);
this->grad(b).multiply(bfs[b-1]->dNdu, Jac);
return true;
}
@ -109,11 +107,11 @@ bool MxFiniteElement::Jacobian (Matrix& Jac, const Matrix& Xnod,
bool MxFiniteElement::Jacobian (Matrix& Jac, Vec3& n,
const Matrix& Xnod,
unsigned short int gBasis,
const std::vector<Matrix>& dNxdu,
const BasisValuesPtrs& bfs,
size_t t1, size_t t2, size_t nBasis,
const Matrix* Xnod2)
{
if (nBasis == 0) nBasis = dNxdu.size();
if (nBasis == 0) nBasis = bfs.size();
const bool separateGeometry = nBasis > this->getNoBasis();
if (separateGeometry) gBasis = nBasis;
@ -123,21 +121,21 @@ bool MxFiniteElement::Jacobian (Matrix& Jac, Vec3& n,
// We are are doing interface terms, evaluate for the second element first
detJxW = utl::Jacobian(Jac,n,this->grad(nBasis+gBasis),
Xnod2 ? *Xnod2 : Xnod,
dNxdu[nBasis+gBasis-1],t1,t2);
bfs[nBasis+gBasis-1]->dNdu,t1,t2);
for (size_t b = 1; b <= nBasis; ++b)
if (b != gBasis)
this->grad(nBasis+b).multiply(dNxdu[nBasis+b-1],Jac);
this->grad(nBasis+b).multiply(bfs[nBasis+b-1]->dNdu,Jac);
}
else if (separateGeometry)
nBasis = this->getNoBasis();
Matrix& dX = separateGeometry ? dummy : this->grad(gBasis);
detJxW = utl::Jacobian(Jac,n,dX,Xnod,dNxdu[gBasis-1],t1,t2);
detJxW = utl::Jacobian(Jac,n,dX,Xnod,bfs[gBasis-1]->dNdu,t1,t2);
for (size_t b = 1; b <= nBasis; ++b)
if (b != gBasis || separateGeometry)
this->grad(b).multiply(dNxdu[b-1],Jac);
this->grad(b).multiply(bfs[b-1]->dNdu,Jac);
return detJxW != 0.0;
}
@ -152,23 +150,21 @@ bool MxFiniteElement::Jacobian (Matrix& Jac, Vec3& n,
bool MxFiniteElement::Hessian (Matrix3D& Hess, const Matrix& Jac,
const Matrix& Xnod,
unsigned short int gBasis,
const std::vector<const BasisFunctionVals*>* bf,
const std::vector<Matrix3D>* d2Nxdu2)
const BasisValuesPtrs& bfs)
{
size_t nBasis = bf ? bf->size() : d2Nxdu2->size();
const bool separateGeometry = nBasis > this->getNoBasis();
const bool separateGeometry = bfs.size() > this->getNoBasis();
bool ok;
if (separateGeometry)
ok = Hess.multiply(Xnod, bf ? bf->back()->d2Ndu2 : d2Nxdu2->back());
ok = Hess.multiply(Xnod, bfs.back()->d2Ndu2);
else
ok = utl::Hessian(Hess, this->hess(gBasis), Jac, Xnod,
bf ? (*bf)[gBasis-1]->d2Ndu2 : (*d2Nxdu2)[gBasis-1],
bfs[gBasis-1]->d2Ndu2,
this->grad(gBasis));
for (size_t b = 1; b <= this->getNoBasis() && ok; b++)
if (b != gBasis || separateGeometry)
ok = utl::Hessian(Hess, this->hess(b), Jac, Xnod,
bf ? (*bf)[b-1]->d2Ndu2 : (*d2Nxdu2)[b-1],
bfs[b-1]->d2Ndu2,
this->grad(b), false);
return ok;

19
src/ASM/FiniteElement.h
View File

@ -14,13 +14,12 @@
#ifndef _FINITE_ELEMENT_H
#define _FINITE_ELEMENT_H
#include "BasisFunctionVals.h"
#include "ItgPoint.h"
#include "MatVec.h"
#include "Vec3.h"
#include "Tensor.h"
struct BasisFunctionVals;
/*!
\brief Class representing a finite element.
@ -121,26 +120,24 @@ public:
//! \param[out] Jac The inverse of the Jacobian matrix
//! \param[in] Xnod Matrix of element nodal coordinates
//! \param[in] gBasis 1-based index of basis representing the geometry
//! \param[in] bf Basis function values and derivatives
//! \param[in] dNxdu First order derivatives of basis functions
//! \param[in] bfs Basis function values and derivatives
bool Jacobian(Matrix& Jac, const Matrix& Xnod,
unsigned short int gBasis,
const std::vector<const BasisFunctionVals*>* bf,
const std::vector<Matrix>* dNxdu = nullptr);
const BasisValuesPtrs& bfs);
//! \brief Sets up the Jacobian matrix of the coordinate mapping on a boundary.
//! \param[out] Jac The inverse of the Jacobian matrix
//! \param[out] n Outward-directed unit normal vector on the boundary
//! \param[in] Xnod Matrix of element nodal coordinates
//! \param[in] gBasis 1-based index of basis representing the geometry
//! \param[in] dNxdu First order derivatives of basis functions
//! \param[in] bfs Derivatives of basis functions
//! \param[in] t1 First parametric tangent direction of the boundary
//! \param[in] t2 Second parametric tangent direction of the boundary
//! \param[in] nBasis Number of basis functions
//! \param[in] Xnod2 Matrix of element nodal coordinates for neighbor element
bool Jacobian(Matrix& Jac, Vec3& n, const Matrix& Xnod,
unsigned short int gBasis,
const std::vector<Matrix>& dNxdu,
const BasisValuesPtrs& bfs,
size_t t1, size_t t2, size_t nBasis = 0,
const Matrix* Xnod2 = nullptr);
@ -149,12 +146,10 @@ public:
//! \param[in] Jac The inverse of the Jacobian matrix
//! \param[in] Xnod Matrix of element nodal coordinates
//! \param[in] gBasis 1-based index of basis representing the geometry
//! \param[in] bf Basis function values and derivatives
//! \param[in] d2Nxdu2 Second order derivatives of basis functions
//! \param[in] bfs Basis function derivatives
bool Hessian(Matrix3D& Hess, const Matrix& Jac, const Matrix& Xnod,
unsigned short int gBasis,
const std::vector<const BasisFunctionVals*>* bf,
const std::vector<Matrix3D>* d2Nxdu2 = nullptr);
const BasisValuesPtrs& bfs);
protected:
//! \brief Returns a reference to the basis function derivatives.

63
src/ASM/LR/ASMu2Dmx.C
View File

@ -426,7 +426,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
fe.u = param[0] = cache.getParam(0,geoEl-1,i);
fe.v = param[1] = cache.getParam(1,geoEl-1,j);
std::vector<const BasisFunctionVals*> bfs(myCache.size());
BasisValuesPtrs bfs(myCache.size());
for (size_t b = 0; b < bfs.size(); ++b) {
bfs[b] = &myCache[b]->getVals(geoEl-1,ig);
if (b < m_basis.size())
@ -435,11 +435,11 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,Xnod,itgBasis,&bfs))
if (!fe.Jacobian(Jac,Xnod,itgBasis,bfs))
continue; // skip singular points
// Compute Hessian of coordinate mapping and 2nd order derivatives
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,itgBasis,&bfs))
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,itgBasis,bfs))
ok = false;
// Compute G-matrix
@ -514,9 +514,8 @@ bool ASMu2Dmx::integrate (Integrand& integrand, int lIndex,
std::map<char,size_t>::const_iterator iit = firstBp.find(lIndex%10);
size_t firstp = iit == firstBp.end() ? 0 : iit->second;
std::vector<Matrix> dNxdu(m_basis.size() + separateGeometry);
BasisValues bfs(m_basis.size() + separateGeometry);
Matrix Xnod, Jac;
Vector Ng;
double param[3] = { 0.0, 0.0, 0.0 };
Vec4 X(param,time.t);
Vec3 normal;
@ -590,25 +589,25 @@ bool ASMu2Dmx::integrate (Integrand& integrand, int lIndex,
if (separateGeometry) {
this->computeBasis(fe.u,fe.v,splinex.back(),els.back()-1,
this->getBasis(ASM::GEOMETRY_BASIS));
SplineUtils::extractBasis(splinex.back(), Ng, dNxdu.back());
SplineUtils::extractBasis(splinex.back(),bfs.back().N,bfs.back().dNdu);
}
// Evaluate basis function derivatives at current integration points
for (size_t b = 0; b < m_basis.size(); ++b) {
this->computeBasis(fe.u, fe.v, splinex[b], els[b]-1, this->getBasis(b+1));
SplineUtils::extractBasis(splinex[b],fe.basis(b+1),dNxdu[b]);
SplineUtils::extractBasis(splinex[b],fe.basis(b+1),bfs[b].dNdu);
}
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,dNxdu,t1,t2))
if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,bfs,t1,t2))
continue; // skip singular points
if (edgeDir < 0)
normal *= -1.0;
// Cartesian coordinates of current integration point
X.assign(Xnod * (separateGeometry ? Ng : fe.basis(itgBasis)));
X.assign(Xnod * (separateGeometry ? bfs.back().N : fe.basis(itgBasis)));
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= dS*wg[i];
@ -777,17 +776,17 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
fe.v = gpar[1][g];
// Evaluate basis function derivatives at current integration points
std::vector<Matrix> dNxdu(2*nB);
BasisValues bfs(2*nB);
for (size_t b = 0; b < nB; b++) {
Go::BasisDerivsSf spline;
this->computeBasis(fe.u+epsu, fe.v+epsv, spline, els[b]-1, m_basis[b].get());
SplineUtils::extractBasis(spline,fe.basis(b+1),dNxdu[b]);
this->computeBasis(fe.u-epsu, fe.v-epsv, spline, els2[b]-1, m_basis[b].get());
SplineUtils::extractBasis(spline,fe.basis(nB+b+1),dNxdu[nB+b]);
this->computeBasis(fe.u+epsu,fe.v+epsv,spline,els[b]-1,m_basis[b].get());
SplineUtils::extractBasis(spline,fe.basis(b+1),bfs[b].dNdu);
this->computeBasis(fe.u-epsu,fe.v-epsv,spline,els2[b]-1,m_basis[b].get());
SplineUtils::extractBasis(spline,fe.basis(nB+b+1),bfs[nB+b].dNdu);
}
// Compute basis function derivatives and the edge normal
if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,dNxdu,t1,t2,nB,&Xnod2))
if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,bfs,t1,t2,nB,&Xnod2))
continue; // skip singular points
if (edgeDir < 0)
@ -834,10 +833,6 @@ bool ASMu2Dmx::evalSolution (Matrix& sField, const Vector& locSol,
if (nPoints != gpar[1].size())
return false;
Vector ptSol;
std::vector<Matrix> dNxdu(m_basis.size()), dNxdX(m_basis.size());
Matrix Jac, Xnod, eSol, ptDer;
Go::BasisPtsSf spline;
std::vector<size_t> nc(nfx.size(), 0);
@ -911,45 +906,47 @@ bool ASMu2Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
this->getElementsAt({gpar[0][i],gpar[1][i]},els,&elem_sizes);
// Evaluate the basis functions at current parametric point
MxFiniteElement fe(elem_sizes,firstIp+i);
std::vector<Matrix> dNxdu(m_basis.size() + separateGeometry);
std::vector<Matrix3D> d2Nxdu2(use2ndDer ? m_basis.size() + separateGeometry : 0);
MxFiniteElement fe(elem_sizes,firstIp+i);
BasisValues bfs(m_basis.size() + separateGeometry);
Matrix Jac, Xnod;
Vector Ng;
Matrix3D Hess;
if (use2ndDer)
for (size_t b = 0; b < d2Nxdu2.size(); ++b) {
for (size_t b = 0; b < bfs.size(); ++b) {
Go::BasisDerivsSf2 spline;
this->computeBasis(gpar[0][i],gpar[1][i],spline,els[b]-1,
b < m_basis.size() ? m_basis[b].get() : geo);
SplineUtils::extractBasis(spline,b < m_basis.size() ? fe.basis(b+1) : Ng,
dNxdu[b],d2Nxdu2[b]);
SplineUtils::extractBasis(spline,
b < m_basis.size() ? fe.basis(b+1) : bfs.back().N,
bfs[b].dNdu,bfs[b].d2Ndu2);
}
else
for (size_t b = 0; b < dNxdu.size(); ++b) {
for (size_t b = 0; b < bfs.size(); ++b) {
Go::BasisDerivsSf spline;
this->computeBasis(gpar[0][i],gpar[1][i],spline,els[b]-1,
b < m_basis.size() ? m_basis[b].get() : geo);
SplineUtils::extractBasis(spline,b < m_basis.size() ? fe.basis(b+1) : Ng,
dNxdu[b]);
SplineUtils::extractBasis(spline,
b < m_basis.size() ? fe.basis(b+1) : bfs.back().N,
bfs[b].dNdu);
}
// Set up control point (nodal) coordinates for current element
if (!this->getElementCoordinates(Xnod,els[itgBasis-1])) return false;
if (!this->getElementCoordinates(Xnod,els[itgBasis-1]))
return false;
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,Xnod,itgBasis,nullptr,&dNxdu))
if (!fe.Jacobian(Jac,Xnod,itgBasis,bfs))
continue; // skip singular points
// Compute Hessian of coordinate mapping and 2nd order derivatives
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,itgBasis,nullptr,&d2Nxdu2))
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,itgBasis,bfs))
return false;
// Cartesian coordinates of current integration point
fe.u = gpar[0][i];
fe.v = gpar[1][i];
utl::Point X4(Xnod * (separateGeometry ? Ng : fe.basis(itgBasis)),{fe.u,fe.v});
utl::Point X4(Xnod * (separateGeometry ? bfs.back().N : fe.basis(itgBasis)),
{fe.u,fe.v});
// Now evaluate the solution field
Vector solPt;

50
src/ASM/LR/ASMu3Dmx.C
View File

@ -451,7 +451,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
fe.v = param[1] = cache.getParam(1,iEl,j);
fe.w = param[2] = cache.getParam(2,iEl,k);
std::vector<const BasisFunctionVals*> bfs(myCache.size());
BasisValuesPtrs bfs(myCache.size());
for (size_t b = 0; b < myCache.size(); ++b) {
bfs[b] = &myCache[b]->getVals(iEl,ig);
if (b < m_basis.size())
@ -460,11 +460,11 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,Xnod,itgBasis,&bfs))
if (!fe.Jacobian(Jac,Xnod,itgBasis,bfs))
continue; // skip singular points
// Compute Hessian of coordinate mapping and 2nd order derivatives
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,itgBasis,&bfs))
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,itgBasis,bfs))
ok = false;
// Compute G-matrix
@ -579,12 +579,11 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex,
fe.v = gpar[1](1);
fe.w = gpar[2](1);
std::vector<Matrix> dNxdu(m_basis.size() + separateGeometry);
BasisValues bfs(m_basis.size() + separateGeometry);
Matrix Xnod, Jac;
double param[3] = { fe.u, fe.v, fe.w };
Vec4 X(param,time.t);
Vector Ng;
Vec3 normal;
double dXidu[3];
@ -645,14 +644,14 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex,
if (separateGeometry)
this->evaluateBasis(els.back()-1, fe.u, fe.v, fe.w,
Ng, dNxdu.back(), ASM::GEOMETRY_BASIS);
bfs.back().N, bfs.back().dNdu, ASM::GEOMETRY_BASIS);
// Fetch basis function derivatives at current integration point
for (size_t b = 1; b <= m_basis.size(); ++b)
this->evaluateBasis(els[b-1]-1, fe, dNxdu[b-1], b);
this->evaluateBasis(els[b-1]-1, fe, bfs[b-1].dNdu, b);
// Compute basis function derivatives and the face normal
if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,dNxdu,t1,t2))
if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,bfs,t1,t2))
continue; // skip singular points
if (faceDir < 0) normal *= -1.0;
@ -662,7 +661,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex,
utl::getGmat(Jac,dXidu,fe.G);
// Cartesian coordinates of current integration point
X.assign(Xnod * (separateGeometry ? Ng : fe.basis(itgBasis)));
X.assign(Xnod * (separateGeometry ? bfs.back().N : fe.basis(itgBasis)));
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= dA*wg[i]*wg[j];
@ -694,9 +693,6 @@ bool ASMu3Dmx::evalSolution (Matrix& sField, const Vector& locSol,
if (nPoints != gpar[1].size() || nPoints != gpar[2].size())
return false;
Vector ptSol;
Matrix Jac, Xnod, eSol, ptDer;
Go::BasisPts spline;
std::vector<size_t> nc(nfx.size(), 0);
@ -770,46 +766,48 @@ bool ASMu3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
this->getElementsAt({gpar[0][i],gpar[1][i],gpar[2][i]},els,&elem_sizes);
// Evaluate the basis functions at current parametric point
MxFiniteElement fe(elem_sizes,firstIp+i);
std::vector<Matrix> dNxdu(m_basis.size() + separateGeometry);
std::vector<Matrix3D> d2Nxdu2(use2ndDer ? m_basis.size() + separateGeometry : 0);
MxFiniteElement fe(elem_sizes,firstIp+i);
BasisValues bfs(m_basis.size() + separateGeometry);
Matrix Jac, Xnod;
Vector Ng;
Matrix3D Hess;
if (use2ndDer)
for (size_t b = 0; b < d2Nxdu2.size(); ++b) {
for (size_t b = 0; b < bfs.size(); ++b) {
Go::BasisDerivs2 spline;
const LR::LRSplineVolume* sv = b < m_basis.size() ? m_basis[b].get() : geo;
sv->computeBasis(gpar[0][i],gpar[1][i],gpar[2][i],spline,els[b]-1);
SplineUtils::extractBasis(spline,b < m_basis.size() ? fe.basis(b+1) : Ng,
dNxdu[b],d2Nxdu2[b]);
SplineUtils::extractBasis(spline,
b < m_basis.size() ? fe.basis(b+1) : bfs.back().N,
bfs[b].dNdu,bfs[b].d2Ndu2);
}
else
for (size_t b = 0; b < dNxdu.size(); ++b) {
for (size_t b = 0; b < bfs.size(); ++b) {
Go::BasisDerivs spline;
const LR::LRSplineVolume* sv = b < m_basis.size() ? m_basis[b].get() : geo;
sv->computeBasis(gpar[0][i],gpar[1][i],gpar[2][i],spline,els[b]-1);
SplineUtils::extractBasis(spline,b < m_basis.size() ? fe.basis(b+1) : Ng,
dNxdu[b]);
SplineUtils::extractBasis(spline,
b < m_basis.size() ? fe.basis(b+1) : bfs.back().N,
bfs[b].dNdu);
}
// Set up control point (nodal) coordinates for current element
if (!this->getElementCoordinates(Xnod,els[itgBasis-1])) return false;
if (!this->getElementCoordinates(Xnod,els[itgBasis-1]))
return false;
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates
if (!fe.Jacobian(Jac,Xnod,itgBasis,nullptr,&dNxdu))
if (!fe.Jacobian(Jac,Xnod,itgBasis,bfs))
continue; // skip singular points
// Compute Hessian of coordinate mapping and 2nd order derivatives
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,itgBasis,nullptr,&d2Nxdu2))
if (use2ndDer && !fe.Hessian(Hess,Jac,Xnod,itgBasis,bfs))
return false;
// Cartesian coordinates of current integration point
fe.u = gpar[0][i];
fe.v = gpar[1][i];
fe.w = gpar[2][i];
utl::Point X4(Xnod * (separateGeometry ? Ng : fe.basis(itgBasis)), {fe.u, fe.v, fe.w});
utl::Point X4(Xnod * (separateGeometry ? bfs.back().N : fe.basis(itgBasis)),
{fe.u, fe.v, fe.w});
// Now evaluate the solution field
Vector solPt;