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

added: FiniteElement::Jacobian for use on boundaries

Browse Source
This commit is contained in:
Arne Morten Kvarving 2023-08-31 09:38:19 +02:00
parent 89c065edc4
commit 95afb36eb2
9 changed files with 170 additions and 166 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

47
Apps/Common/Test/TestCompatibleOperators.C
View File

@ -22,29 +22,30 @@ const auto check_matrix_equal = [](const Matrix& A, const DoubleVec& B)
ASSERT_NEAR(A(i,j), B[i-1][j-1], 1e-13); ASSERT_NEAR(A(i,j), B[i-1][j-1], 1e-13);
}; };
static MxFiniteElement getFE() class MyFiniteElement : public MxFiniteElement
{ {
MxFiniteElement fe({6,6,4}); public:
for (size_t i = 1; i <= 2; ++i) { MyFiniteElement() : MxFiniteElement({6,6,4})
fe.grad(i).resize(6,2); {
fe.basis(i).resize(6); for (size_t i = 1; i <= 2; ++i) {
for (size_t j = 1; j <= 6; ++j) { this->grad(i).resize(6,2);
fe.grad(i)(j,1) = 1.0+12*(i-1)+2*(j-1); this->basis(i).resize(6);
fe.grad(i)(j,2) = 2.0+12*(i-1)+2*(j-1); for (size_t j = 1; j <= 6; ++j) {
fe.basis(i)(j) = j + 6*(i-1); this->grad(i)(j,1) = 1.0+12*(i-1)+2*(j-1);
this->grad(i)(j,2) = 2.0+12*(i-1)+2*(j-1);
this->basis(i)(j) = j + 6*(i-1);
}
}
this->grad(3).resize(4,2);
this->basis(3).resize(4);
for (size_t j = 1; j <= 4; ++j) {
this->grad(3)(j,1) = 1.0 + 14.0 + 2*(j-1);
this->grad(3)(j,2) = 2.0 + 14.0 + 2*(j-1);
this->basis(3)(j) = 12 + j;
} }
} }
};
fe.grad(3).resize(4,2);
fe.basis(3).resize(4);
for (size_t j = 1; j <= 4; ++j) {
fe.grad(3)(j,1) = 1.0 + 14.0 + 2*(j-1);
fe.grad(3)(j,2) = 2.0 + 14.0 + 2*(j-1);
fe.basis(3)(j) = 12 + j;
}
return fe;
}
/* /*
@ -118,7 +119,7 @@ TESTI(TestEqualOrderOperators, Gradient)
TEST(TestCompatibleOperators, Laplacian) TEST(TestCompatibleOperators, Laplacian)
{ {
MxFiniteElement fe = getFE(); MyFiniteElement fe;
std::vector<Matrix> EM(3); std::vector<Matrix> EM(3);
EM[1].resize(6,6); EM[1].resize(6,6);
@ -174,7 +175,7 @@ TEST(TestCompatibleOperators, Laplacian)
TEST(TestCompatibleOperators, Mass) TEST(TestCompatibleOperators, Mass)
{ {
MxFiniteElement fe = getFE(); MyFiniteElement fe;
std::vector<Matrix> EM_vec(3); std::vector<Matrix> EM_vec(3);
EM_vec[1].resize(6,6); EM_vec[1].resize(6,6);
@ -204,7 +205,7 @@ TEST(TestCompatibleOperators, Mass)
TEST(TestCompatibleOperators, Source) TEST(TestCompatibleOperators, Source)
{ {
MxFiniteElement fe = getFE(); MyFiniteElement fe;
Vectors EV_scalar(3); Vectors EV_scalar(3);
EV_scalar[1].resize(6); EV_scalar[1].resize(6);

69
src/ASM/ASMs2Dmx.C
View File

@ -691,41 +691,36 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
for (int i = 0; i < nGauss && ok; i++, ip++, fe.iGP++) for (int i = 0; i < nGauss && ok; i++, ip++, fe.iGP++)
{ {
// Local element coordinates and parameter values // Local element coordinates and parameter values
// of current integration point // of current integration point
if (gpar[0].size() > 1) if (gpar[0].size() > 1)
{ {
fe.xi = xg[i]; fe.xi = xg[i];
fe.u = param[0] = gpar[0](i+1,i1-p1+1); fe.u = param[0] = gpar[0](i+1,i1-p1+1);
} }
if (gpar[1].size() > 1) if (gpar[1].size() > 1)
{ {
fe.eta = xg[i]; fe.eta = xg[i];
fe.v = param[1] = gpar[1](i+1,i2-p2+1); fe.v = param[1] = gpar[1](i+1,i2-p2+1);
} }
// Fetch basis function derivatives at current integration point // Fetch basis function derivatives at current integration point
for (size_t b = 0; b < m_basis.size(); ++b) 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),dNxdu[b]);
// Compute Jacobian inverse of the coordinate mapping and // Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates // basis function derivatives w.r.t. Cartesian coordinates
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(itgBasis),Xnod, if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,dNxdu,t1,t2))
dNxdu[itgBasis-1],t1,t2); continue; // skip singular points
if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 0; b < m_basis.size(); ++b) if (edgeDir < 0) normal *= -1.0;
if (b != (size_t)itgBasis-1)
fe.grad(b+1).multiply(dNxdu[b],Jac);
if (edgeDir < 0) normal *= -1.0; // Cartesian coordinates of current integration point
X.assign(Xnod * fe.basis(itgBasis));
// Cartesian coordinates of current integration point // Evaluate the integrand and accumulate element contributions
X.assign(Xnod * fe.basis(itgBasis)); fe.detJxW *= dS*wg[i];
ok = integrand.evalBouMx(*A,fe,time,X,normal);
// Evaluate the integrand and accumulate element contributions
fe.detJxW *= dS*wg[i];
ok = integrand.evalBouMx(*A,fe,time,X,normal);
} }
// Finalize the element quantities // Finalize the element quantities
@ -777,9 +772,10 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
const int p2 = surf->order_v(); const int p2 = surf->order_v();
const int n1 = surf->numCoefs_u(); const int n1 = surf->numCoefs_u();
const int n2 = surf->numCoefs_v(); const int n2 = surf->numCoefs_v();
const size_t nB = m_basis.size();
std::vector<size_t> elem_sizes2(elem_size); std::vector<size_t> elem_sizes2(elem_size);
std::copy(elem_size.begin(), elem_size.end(), std::back_inserter(elem_sizes2)); std::copy(elem_size.begin(),elem_size.end(),std::back_inserter(elem_sizes2));
MxFiniteElement fe(elem_sizes2); MxFiniteElement fe(elem_sizes2);
Matrix dNdu, Xnod, Jac; Matrix dNdu, Xnod, Jac;
@ -877,27 +873,18 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
} }
// Fetch basis function derivatives at current integration point // Fetch basis function derivatives at current integration point
Matrices dNxdu(m_basis.size()*2); Matrices dNxdu(2*nB);
for (size_t b = 0; b < m_basis.size(); ++b) { for (size_t b = 1; b <= nB; b++) {
Go::BasisDerivsSf spline; Go::BasisDerivsSf spline;
this->getBasis(b+1)->computeBasis(fe.u, fe.v, spline, edgeDir < 0); this->getBasis(b)->computeBasis(fe.u, fe.v, spline, edgeDir < 0);
SplineUtils::extractBasis(spline, fe.basis(b+1), dNxdu[b]); SplineUtils::extractBasis(spline, fe.basis(b), dNxdu[b-1]);
this->getBasis(b+1)->computeBasis(fe.u, fe.v, spline, edgeDir > 0); this->getBasis(b)->computeBasis(fe.u, fe.v, spline, edgeDir > 0);
SplineUtils::extractBasis(spline, fe.basis(b+1+m_basis.size()), SplineUtils::extractBasis(spline, fe.basis(nB+b), dNxdu[nB+b-1]);
dNxdu[b+m_basis.size()]);
} }
// Compute basis function derivatives and the edge normal // Compute basis function derivatives and the edge normal
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(itgBasis+m_basis.size()),Xnod, if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,dNxdu,t1,t2,nB))
dNxdu[itgBasis-1+m_basis.size()],t1,t2); continue; // skip singular points
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(itgBasis),Xnod,
dNxdu[itgBasis-1],t1,t2);
if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 0; b < m_basis.size(); ++b)
if (b != (size_t)itgBasis-1) {
fe.grad(b+1).multiply(dNxdu[b],Jac);
fe.grad(b+1+m_basis.size()).multiply(dNxdu[b+m_basis.size()],Jac);
}
if (edgeDir < 0) normal *= -1.0; if (edgeDir < 0) normal *= -1.0;

35
src/ASM/ASMs2DmxLag.C
View File

@ -415,12 +415,12 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, int lIndex,
for (int i = 0; i < nGauss && ok; i++, fe.iGP++) for (int i = 0; i < nGauss && ok; i++, fe.iGP++)
{ {
// Gauss point coordinates along the edge // Gauss point coordinates along the edge
xi[t1-1] = edgeDir < 0 ? -1.0 : 1.0; xi[t1-1] = edgeDir < 0 ? -1.0 : 1.0;
xi[t2-1] = xg[i]; xi[t2-1] = xg[i];
// Compute the basis functions and their derivatives, using // Compute the basis functions and their derivatives, using
// tensor product of one-dimensional Lagrange polynomials // tensor product of one-dimensional Lagrange polynomials
for (size_t b = 0; b < nxx.size(); ++b) 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),dNxdu[b],
elem_sizes[b][0],xi[0], elem_sizes[b][0],xi[0],
@ -428,23 +428,18 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, int lIndex,
ok = false; ok = false;
// Compute basis function derivatives and the edge normal // Compute basis function derivatives and the edge normal
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(itgBasis),Xnod, if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,dNxdu,t1,t2))
dNxdu[itgBasis-1],t1,t2); continue; // skip singular points
if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 0; b < nxx.size(); ++b) if (edgeDir < 0) normal *= -1.0;
if (b != (size_t)itgBasis-1)
fe.grad(b+1).multiply(dNxdu[b],Jac);
if (edgeDir < 0) normal *= -1.0; // Cartesian coordinates of current integration point
X.assign(Xnod * fe.basis(itgBasis));
// Cartesian coordinates of current integration point // Evaluate the integrand and accumulate element contributions
X.assign(Xnod * fe.basis(itgBasis)); fe.detJxW *= wg[i];
if (ok && !integrand.evalBouMx(*A,fe,time,X,normal))
// Evaluate the integrand and accumulate element contributions ok = false;
fe.detJxW *= wg[i];
if (ok && !integrand.evalBouMx(*A,fe,time,X,normal))
ok = false;
} }
// Finalize the element quantities // Finalize the element quantities
@ -453,7 +448,7 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, int lIndex,
// Assembly of global system integral // Assembly of global system integral
if (ok && !glInt.assemble(A->ref(),fe.iel)) if (ok && !glInt.assemble(A->ref(),fe.iel))
ok = false; ok = false;
A->destruct(); A->destruct();

35
src/ASM/ASMs3Dmx.C
View File

@ -822,13 +822,8 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
// Compute Jacobian inverse of the coordinate mapping and // Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates // basis function derivatives w.r.t. Cartesian coordinates
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(itgBasis),Xnod, if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,dNxdu,t1,t2))
dNxdu[itgBasis-1],t1,t2); continue; // skip singular points
if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 0; b < m_basis.size(); ++b)
if (b != (size_t)itgBasis-1)
fe.grad(b+1).multiply(dNxdu[b],Jac);
if (faceDir < 0) normal *= -1.0; if (faceDir < 0) normal *= -1.0;
@ -894,6 +889,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
const int n3 = svol->numCoefs(2); const int n3 = svol->numCoefs(2);
const int nel1 = n1 - p1 + 1; const int nel1 = n1 - p1 + 1;
const int nel2 = n2 - p2 + 1; const int nel2 = n2 - p2 + 1;
const size_t nB = m_basis.size();
std::vector<size_t> elem_sizes2(elem_size); std::vector<size_t> elem_sizes2(elem_size);
std::copy(elem_size.begin(), elem_size.end(), std::back_inserter(elem_sizes2)); std::copy(elem_size.begin(), elem_size.end(), std::back_inserter(elem_sizes2));
@ -1028,27 +1024,18 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
} }
// Fetch basis function derivatives at current integration point // Fetch basis function derivatives at current integration point
Matrices dNxdu(m_basis.size()*2); Matrices dNxdu(2*nB);
for (size_t b = 0; b < m_basis.size(); ++b) { for (size_t b = 1; b <= nB; b++) {
Go::BasisDerivs spline; Go::BasisDerivs spline;
this->getBasis(b+1)->computeBasis(fe.u, fe.v, fe.w, spline, faceDir < 0); this->getBasis(b)->computeBasis(fe.u, fe.v, fe.w, spline, faceDir < 0);
SplineUtils::extractBasis(spline, fe.basis(b+1), dNxdu[b]); SplineUtils::extractBasis(spline, fe.basis(b), dNxdu[b-1]);
this->getBasis(b+1)->computeBasis(fe.u, fe.v, fe.w, spline, faceDir > 0); this->getBasis(b)->computeBasis(fe.u, fe.v, fe.w, spline, faceDir > 0);
SplineUtils::extractBasis(spline, fe.basis(b+1+m_basis.size()), SplineUtils::extractBasis(spline, fe.basis(nB+b), dNxdu[nB+b-1]);
dNxdu[b+m_basis.size()]);
} }
// Compute basis function derivatives and the edge normal // Compute basis function derivatives and the edge normal
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(itgBasis+m_basis.size()), if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,dNxdu,t1,t2,nB))
Xnod,dNxdu[itgBasis-1+m_basis.size()],t1,t2); continue; // skip singular points
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(itgBasis),Xnod,
dNxdu[itgBasis-1],t1,t2);
if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 0; b < m_basis.size(); ++b)
if (b != (size_t)itgBasis-1) {
fe.grad(b+1).multiply(dNxdu[b],Jac);
fe.grad(b+1+m_basis.size()).multiply(dNxdu[b+m_basis.size()],Jac);
}
if (faceDir < 0) normal *= -1.0; if (faceDir < 0) normal *= -1.0;

15
src/ASM/ASMs3DmxLag.C
View File

@ -493,19 +493,14 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, int lIndex,
elem_sizes[b][2],xi[2])) elem_sizes[b][2],xi[2]))
ok = false; ok = false;
// Compute basis function derivatives and the edge normal // Compute basis function derivatives and the edge normal
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(itgBasis),Xnod, if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,dNxdu,t1,t2))
dNxdu[itgBasis-1],t1,t2); continue; // skip singular points
if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 0; b < nxx.size(); ++b) if (faceDir < 0) normal *= -1.0;
if (b != (size_t)itgBasis-1)
fe.grad(b+1).multiply(dNxdu[b],Jac);
if (faceDir < 0) normal *= -1.0;
// Cartesian coordinates of current integration point // Cartesian coordinates of current integration point
X.assign(Xnod * fe.basis(itgBasis)); X.assign(Xnod * fe.basis(itgBasis));
// Evaluate the integrand and accumulate element contributions // Evaluate the integrand and accumulate element contributions
fe.detJxW *= wg[i]*wg[j]; fe.detJxW *= wg[i]*wg[j];

50
src/ASM/FiniteElement.C
View File

@ -98,6 +98,51 @@ bool MxFiniteElement::Jacobian (Matrix& Jac, const Matrix& Xnod,
} }
/*!
If the input argument \a nBasis is half (or less than half) of the size of
the internal basis function value array, it is used to flag that we are doing
element interface terms, and the basis function values of both elements
sharing the interface are stored internally. If \a nBasis is zero (default),
it is reset to the size of the \a dNxdu argument, which will fit normal cases.
*/
bool MxFiniteElement::Jacobian (Matrix& Jac, Vec3& n,
const Matrix& Xnod,
unsigned short int gBasis,
const std::vector<Matrix>& dNxdu,
size_t t1, size_t t2, size_t nBasis,
const Matrix* Xnod2)
{
if (nBasis == 0) nBasis = dNxdu.size();
const bool separateGeometry = nBasis > this->getNoBasis();
if (separateGeometry) gBasis = nBasis;
Matrix dummy;
if (2*nBasis <= this->getNoBasis())
{
// 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);
for (size_t b = 1; b <= nBasis; ++b)
if (b != gBasis)
this->grad(nBasis+b).multiply(dNxdu[nBasis+b-1],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);
for (size_t b = 1; b <= nBasis; ++b)
if (b != gBasis || separateGeometry)
this->grad(b).multiply(dNxdu[b-1],Jac);
return detJxW != 0.0;
}
/*! /*!
This method also calculates the second-derivatives of the basis functions This method also calculates the second-derivatives of the basis functions
with respect to the Cartesian coordinates, using the same geometry mapping with respect to the Cartesian coordinates, using the same geometry mapping
@ -114,12 +159,11 @@ bool MxFiniteElement::Hessian (Matrix3D& Hess, const Matrix& Jac,
const bool separateGeometry = nBasis > this->getNoBasis(); const bool separateGeometry = nBasis > this->getNoBasis();
bool ok; bool ok;
if (separateGeometry) if (separateGeometry)
ok = Hess.multiply(Xnod, (bf ? bf->back()->d2Ndu2 : d2Nxdu2->back())); ok = Hess.multiply(Xnod, bf ? bf->back()->d2Ndu2 : d2Nxdu2->back());
else { else
ok = utl::Hessian(Hess, this->hess(gBasis), Jac, Xnod, ok = utl::Hessian(Hess, this->hess(gBasis), Jac, Xnod,
bf ? (*bf)[gBasis-1]->d2Ndu2 : (*d2Nxdu2)[gBasis-1], bf ? (*bf)[gBasis-1]->d2Ndu2 : (*d2Nxdu2)[gBasis-1],
this->grad(gBasis)); this->grad(gBasis));
}
for (size_t b = 1; b <= this->getNoBasis() && ok; b++) for (size_t b = 1; b <= this->getNoBasis() && ok; b++)
if (b != gBasis || separateGeometry) if (b != gBasis || separateGeometry)

36
src/ASM/FiniteElement.h
View File

@ -46,17 +46,17 @@ public:
//! \brief Returns a const reference to the basis function derivatives. //! \brief Returns a const reference to the basis function derivatives.
virtual const Matrix& grad(char) const { return dNdX; } virtual const Matrix& grad(char) const { return dNdX; }
//! \brief Returns a reference to the basis function derivatives.
virtual Matrix& grad(char) { return dNdX; }
//! \brief Returns a const reference to the basis function 2nd-derivatives. //! \brief Returns a const reference to the basis function 2nd-derivatives.
virtual const Matrix3D& hess(char) const { return d2NdX2; } virtual const Matrix3D& hess(char) const { return d2NdX2; }
//! \brief Returns a reference to the basis function 2nd-derivatives.
virtual Matrix3D& hess(char) { return d2NdX2; }
//! \brief Returns a const reference to the basis function 3nd-derivatives. //! \brief Returns a const reference to the basis function 3nd-derivatives.
virtual const Matrix4D& hess2(char) const { return d3NdX3; } virtual const Matrix4D& hess2(char) const { return d3NdX3; }
protected: protected:
//! \brief Returns a reference to the basis function derivatives.
virtual Matrix& grad(char) { return dNdX; }
//! \brief Returns a reference to the basis function 2nd-derivatives.
virtual Matrix3D& hess(char) { return d2NdX2; }
//! \brief Writes the finite element object to the given output stream. //! \brief Writes the finite element object to the given output stream.
virtual std::ostream& write(std::ostream& os) const; virtual std::ostream& write(std::ostream& os) const;
@ -112,13 +112,8 @@ public:
//! \brief Returns a const reference to the basis function derivatives. //! \brief Returns a const reference to the basis function derivatives.
virtual const Matrix& grad(char b) const { return b < 2 ? dNdX : dMdX[b-2]; } virtual const Matrix& grad(char b) const { return b < 2 ? dNdX : dMdX[b-2]; }
//! \brief Returns a reference to the basis function derivatives.
virtual Matrix& grad(char b) { return b < 2 ? dNdX : dMdX[b-2]; }
//! \brief Returns a const reference to the basis function 2nd-derivatives. //! \brief Returns a const reference to the basis function 2nd-derivatives.
virtual const Matrix3D& hess(char b) const { return b < 2 ? d2NdX2 : d2MdX2[b-2]; } virtual const Matrix3D& hess(char b) const { return b < 2 ? d2NdX2 : d2MdX2[b-2]; }
//! \brief Returns a reference to the basis function 2nd-derivatives.
virtual Matrix3D& hess(char b) { return b < 2 ? d2NdX2 : d2MdX2[b-2]; }
//! \brief Returns a const reference to the basis function 3rd-derivatives. //! \brief Returns a const reference to the basis function 3rd-derivatives.
virtual const Matrix4D& hess2(char b) const { return b < 2 ? d3NdX3 : d3MdX3[b-2]; } virtual const Matrix4D& hess2(char b) const { return b < 2 ? d3NdX3 : d3MdX3[b-2]; }
@ -133,6 +128,22 @@ public:
const std::vector<const BasisFunctionVals*>* bf, const std::vector<const BasisFunctionVals*>* bf,
const std::vector<Matrix>* dNxdu = nullptr); const std::vector<Matrix>* dNxdu = nullptr);
//! \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] 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,
size_t t1, size_t t2, size_t nBasis = 0,
const Matrix* Xnod2 = nullptr);
//! \brief Sets up the Hessian matrix of the coordinate mapping. //! \brief Sets up the Hessian matrix of the coordinate mapping.
//! \param[out] Hess The Hessian matrix //! \param[out] Hess The Hessian matrix
//! \param[in] Jac The inverse of the Jacobian matrix //! \param[in] Jac The inverse of the Jacobian matrix
@ -146,6 +157,11 @@ public:
const std::vector<Matrix3D>* d2Nxdu2 = nullptr); const std::vector<Matrix3D>* d2Nxdu2 = nullptr);
protected: protected:
//! \brief Returns a reference to the basis function derivatives.
virtual Matrix& grad(char b) { return b < 2 ? dNdX : dMdX[b-2]; }
//! \brief Returns a reference to the basis function 2nd-derivatives.
virtual Matrix3D& hess(char b) { return b < 2 ? d2NdX2 : d2MdX2[b-2]; }
//! \brief Writes the finite element object to the given output stream. //! \brief Writes the finite element object to the given output stream.
virtual std::ostream& write(std::ostream& os) const; virtual std::ostream& write(std::ostream& os) const;

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

@ -574,13 +574,8 @@ bool ASMu2Dmx::integrate (Integrand& integrand, int lIndex,
// Compute Jacobian inverse of the coordinate mapping and // Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates // basis function derivatives w.r.t. Cartesian coordinates
fe.detJxW = utl::Jacobian(Jac,normal,fe.grad(itgBasis),Xnod, if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,dNxdu,t1,t2))
dNxdu[itgBasis-1],t1,t2); continue; // skip singular points
if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 1; b <= m_basis.size(); ++b)
if ((int)b != itgBasis)
fe.grad(b).multiply(dNxdu[b-1],Jac);
if (edgeDir < 0) if (edgeDir < 0)
normal *= -1.0; normal *= -1.0;
@ -632,6 +627,8 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
const double* wg = GaussQuadrature::getWeight(nGP); const double* wg = GaussQuadrature::getWeight(nGP);
if (!xg || !wg) return false; if (!xg || !wg) return false;
const size_t nB = m_basis.size();
Matrix Xnod, Jac; Matrix Xnod, Jac;
Vec4 X(nullptr,time.t); Vec4 X(nullptr,time.t);
Vec3 normal; Vec3 normal;
@ -747,31 +744,18 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
fe.v = gpar[1][g]; fe.v = gpar[1][g];
// Evaluate basis function derivatives at current integration points // Evaluate basis function derivatives at current integration points
std::vector<Matrix> dNxdu(m_basis.size()*2); std::vector<Matrix> dNxdu(2*nB);
for (size_t b=0; b < m_basis.size(); ++b) { for (size_t b = 0; b < nB; b++) {
Go::BasisDerivsSf spline; Go::BasisDerivsSf spline;
this->computeBasis(fe.u+epsu, fe.v+epsv, spline, els[b]-1, m_basis[b].get()); 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]); 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()); this->computeBasis(fe.u-epsu, fe.v-epsv, spline, els2[b]-1, m_basis[b].get());
SplineUtils::extractBasis(spline,fe.basis(b+1+m_basis.size()), SplineUtils::extractBasis(spline,fe.basis(nB+b+1),dNxdu[nB+b]);
dNxdu[b+m_basis.size()]);
} }
// Compute Jacobian inverse of the coordinate mapping and // Compute basis function derivatives and the edge normal
// basis function derivatives w.r.t. Cartesian coordinates if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,dNxdu,t1,t2,nB,&Xnod2))
fe.detJxW = utl::Jacobian(Jac2,normal, continue; // skip singular points
fe.grad(itgBasis+m_basis.size()),Xnod2,
dNxdu[itgBasis-1+m_basis.size()],t1,t2);
fe.detJxW = utl::Jacobian(Jac,normal,
fe.grad(itgBasis),Xnod,
dNxdu[itgBasis-1],t1,t2);
if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 1; b <= m_basis.size(); ++b)
if ((int)b != itgBasis) {
fe.grad(b).multiply(dNxdu[b-1],Jac);
fe.grad(b+m_basis.size()).multiply(dNxdu[b-1+m_basis.size()],Jac);
}
if (edgeDir < 0) if (edgeDir < 0)
normal *= -1.0; normal *= -1.0;

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

@ -376,6 +376,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
{ {
// --- Compute average value of basis functions over the element ----- // --- Compute average value of basis functions over the element -----
Matrix dNdX;
fe.Navg.resize(elem_sizes[0],true); fe.Navg.resize(elem_sizes[0],true);
double vol = 0.0; double vol = 0.0;
size_t jp = 0; size_t jp = 0;
@ -388,8 +389,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
// Compute Jacobian determinant of coordinate mapping // Compute Jacobian determinant of coordinate mapping
// and multiply by weight of current integration point // and multiply by weight of current integration point
double detJac = utl::Jacobian(Jac,fe.grad(itgBasis), double detJac = utl::Jacobian(Jac,dNdX,Xnod,bfs.dNdu,false);
Xnod,bfs.dNdu,false);
double weight = dV*wg[0][i]*wg[1][j]*wg[2][k]; double weight = dV*wg[0][i]*wg[1][j]*wg[2][k];
// Numerical quadrature // Numerical quadrature
@ -634,13 +634,8 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex,
this->evaluateBasis(iEl, fe, dNxdu[b-1], b); this->evaluateBasis(iEl, fe, dNxdu[b-1], b);
// Compute basis function derivatives and the face normal // Compute basis function derivatives and the face normal
fe.detJxW = utl::Jacobian(Jac, normal, fe.grad(itgBasis), Xnod, if (!fe.Jacobian(Jac,normal,Xnod,itgBasis,dNxdu,t1,t2))
dNxdu[itgBasis-1], t1, t2); continue; // skip singular points
if (fe.detJxW == 0.0) continue; // skip singular points
for (size_t b = 1; b <= m_basis.size(); ++b)
if ((int)b != itgBasis)
fe.grad(b).multiply(dNxdu[b-1],Jac);
if (faceDir < 0) normal *= -1.0; if (faceDir < 0) normal *= -1.0;