MxFiniteElement::Jacobian

support for a separate geometry basis. checks length
of the passed vector, if greater than number of bases,
use the last entry for the geometry mapping
This commit is contained in:
Arne Morten Kvarving 2023-08-31 09:34:19 +02:00
parent 004d1927ae
commit 28d547c3e4

View File

@ -80,18 +80,19 @@ bool MxFiniteElement::Jacobian (Matrix& Jac, const Matrix& Xnod,
const std::vector<const BasisFunctionVals*>* bf,
const std::vector<Matrix>* dNxdu)
{
detJxW = utl::Jacobian(Jac,
gBasis > 1 ? dMdX[gBasis-2] : dNdX, Xnod,
bf ? (*bf)[gBasis-1]->dNdu : (*dNxdu)[gBasis-1]);
size_t nBasis = bf ? bf->size() : dNxdu->size();
const bool separateGeometry = nBasis > this->getNoBasis();
if (separateGeometry)
gBasis = nBasis;
detJxW = utl::Jacobian(Jac, this->grad(gBasis), Xnod,
bf ? (*bf)[gBasis-1]->dNdu : (*dNxdu)[gBasis-1],
!separateGeometry);
if (detJxW == 0.0) return false; // singular point
if (gBasis > 1)
dNdX.multiply(bf ? bf->front()->dNdu : dNxdu->front(),Jac);
size_t nBasis = bf ? bf->size() : dNxdu->size();
for (size_t b = 2; b <= nBasis; b++)
if (b != gBasis)
dMdX[b-2].multiply(bf ? (*bf)[b-1]->dNdu : (*dNxdu)[b-1], Jac);
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);
return true;
}