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

Added: Class ItgPoint representing an integration/evaluation point,

Browse Source 
the class FiniteElement then inherits ItgPoint.
Changed: All valueFE methods use ItgPoint instead of FiniteElement.
This commit is contained in:
Knut Morten Okstad 2019-10-02 19:33:15 +02:00
parent a02a9a0251
commit 70176f1408
43 changed files with 557 additions and 540 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists.txt
View File

@ -165,7 +165,7 @@ else()
src/ASM/Integrand.h src/ASM/Lagrange.h
src/ASM/LocalIntegral.h src/ASM/SAMpatch.h
src/ASM/TimeDomain.h src/ASM/ASMs?D.h src/ASM/ASM?D.h
src/ASM/DomainDecomposition.h
src/ASM/DomainDecomposition.h src/ASM/ItgPoint.h
src/LinAlg/*.h src/SIM/*.h src/Utility/*.h
3rdparty/*.h
${CMAKE_BINARY_DIR}/IFEM.h)

8
src/ASM/ASMs1D.C
View File

@ -1560,14 +1560,12 @@ bool ASMs1D::evalProjSolution (Matrix& sField, const Vector& locSol,
// Evaluate the projected solution field at each point
Vector vals;
FiniteElement fe;
sField.resize(f->getNoFields(),gpar.size());
size_t ipt = 0;
for (size_t i = 0; i < gpar.size(); i++)
for (double u : gpar)
{
fe.u = gpar[i];
f->valueFE(fe,vals);
f->valueFE(ItgPoint(u),vals);
sField.fillColumn(++ipt,vals);
}
@ -1861,7 +1859,7 @@ bool ASMs1D::evaluate (const FunctionBase* func, RealArray& values,
Fields* ASMs1D::getProjectedFields (const Vector& coefs, size_t) const
{
if (proj == curv)
if (proj == curv || this->getNoProjectionNodes() == 0)
return nullptr;
size_t ncmp = coefs.size() / this->getNoProjectionNodes();

13
src/ASM/ASMs2D.C
View File

@ -2753,16 +2753,13 @@ bool ASMs2D::evalProjSolution (Matrix& sField, const Vector& locSol,
// Evaluate the projected solution field at each point
Vector vals;
FiniteElement fe;
sField.resize(f->getNoFields(),gpar[0].size()*gpar[1].size());
size_t ipt = 0;
for (size_t j = 0; j < gpar[1].size(); j++)
for (size_t i = 0; i < gpar[0].size(); i++)
for (double v : gpar[1])
for (double u : gpar[0])
{
fe.u = gpar[0][i];
fe.v = gpar[1][j];
f->valueFE(fe,vals);
f->valueFE(ItgPoint(u,v),vals);
sField.fillColumn(++ipt,vals);
}
@ -3126,7 +3123,7 @@ int ASMs2D::getCorner (int I, int J, int basis) const
Fields* ASMs2D::getProjectedFields (const Vector& coefs, size_t) const
{
if (proj == this->getBasis(1))
if (proj == this->getBasis(1) || this->getNoProjectionNodes() == 0)
return nullptr;
size_t ncmp = coefs.size() / this->getNoProjectionNodes();
@ -3142,6 +3139,8 @@ Fields* ASMs2D::getProjectedFields (const Vector& coefs, size_t) const
size_t ASMs2D::getNoProjectionNodes () const
{
if (!proj) return 0;
return proj->numCoefs_u() * proj->numCoefs_v();
}

18
src/ASM/ASMs2Drecovery.C
View File

@ -15,8 +15,7 @@
#include "GoTools/geometry/SurfaceInterpolator.h"
#include "ASMs2D.h"
#include "ASMs2Dmx.h"
#include "FiniteElement.h"
#include "ItgPoint.h"
#include "Field.h"
#include "IntegrandBase.h"
#include "CoordinateMapping.h"
@ -485,16 +484,11 @@ bool ASMs2D::evaluate (const Field* field, RealArray& vec, int basisNum) const
// Evaluate the result field at all sampling points.
// Note: it is here assumed that *basis and *this have spline bases
// defined over the same parameter domain.
Vector sValues(gpar[0].size()*gpar[1].size());
Vector::iterator it=sValues.begin();
for (size_t j=0;j<gpar[1].size();++j) {
FiniteElement fe;
fe.v = gpar[1][j];
for (size_t i=0;i<gpar[0].size();++i) {
fe.u = gpar[0][i];
*it++ = field->valueFE(fe);
}
}
Vector sValues;
sValues.reserve(gpar[0].size()*gpar[1].size());
for (double v : gpar[1])
for (double u : gpar[0])
sValues.push_back(field->valueFE(ItgPoint(u,v)));
Go::SplineSurface* surf = this->getBasis(basisNum);

21
src/ASM/ASMs3Drecovery.C
View File

@ -15,7 +15,7 @@
#include "GoTools/trivariate/VolumeInterpolator.h"
#include "ASMs3D.h"
#include "FiniteElement.h"
#include "ItgPoint.h"
#include "Field.h"
#include "CoordinateMapping.h"
#include "GaussQuadrature.h"
@ -300,19 +300,12 @@ bool ASMs3D::evaluate (const Field* field, RealArray& vec, int basisNum) const
// Evaluate the result field at all sampling points.
// Note: it is here assumed that *basis and *this have spline bases
// defined over the same parameter domain.
Vector sValues(gpar[0].size()*gpar[1].size()*gpar[2].size());
Vector::iterator it=sValues.begin();
for (size_t l=0;l<gpar[2].size();++l) {
FiniteElement fe;
fe.w = gpar[2][l];
for (size_t j=0;j<gpar[1].size();++j) {
fe.v = gpar[1][j];
for (size_t i=0;i<gpar[0].size();++i) {
fe.u = gpar[0][i];
*it++ = field->valueFE(fe);
}
}
}
Vector sValues;
sValues.reserve(gpar[0].size()*gpar[1].size()*gpar[2].size());
for (double w : gpar[2])
for (double v : gpar[1])
for (double u : gpar[0])
sValues.push_back(field->valueFE(ItgPoint(u,v,w)));
Go::SplineVolume* svol = this->getBasis(basisNum);

25
src/ASM/Field.h
View File

@ -18,7 +18,7 @@
#include <string>
class ASMbase;
class FiniteElement;
class ItgPoint;
class Vec4;
@ -54,12 +54,12 @@ public:
//==============================
//! \brief Computes the value in a given node/control point.
//! \param[in] node Node number
//! \param[in] node 1-based node/control point index
virtual double valueNode(size_t node) const = 0;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element quantities
virtual double valueFE(const FiniteElement& fe) const = 0;
//! \param[in] x Local coordinate of evaluation point
virtual double valueFE(const ItgPoint& x) const = 0;
//! \brief Computes the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point
@ -71,21 +71,20 @@ public:
virtual bool valueGrid(RealArray& val, const int* npe) const { return false; }
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element quantities
//! \param[in] x Local coordinate of evaluation point
//! \param[out] grad Gradient of solution in a given local coordinate
virtual bool gradFE(const FiniteElement& fe, Vector& grad) const = 0;
//! \brief Computes the hessian for a given local coordinate.
//! \param[in] fe Finite element quantities
//! \param[out] d2UdX2 Hessian of solution in a given local coordinate
virtual bool hessianFE(const FiniteElement& fe, Matrix& d2UdX2) const
{ return false; }
virtual bool gradFE(const ItgPoint& x, Vector& grad) const = 0;
//! \brief Computes the gradient for a given global/physical coordinate.
//! \param[in] x Global coordinate
//! \param[in] x Global/physical coordinate for point
//! \param[out] grad Gradient of solution in a given global coordinate
virtual bool gradCoor(const Vec4& x, Vector& grad) const { return false; }
//! \brief Computes the hessian for a given local coordinate.
//! \param[in] x Local coordinate of evaluation point
//! \param[out] H Hessian of solution in a given local coordinate
virtual bool hessianFE(const ItgPoint& x, Matrix& H) const { return false; }
protected:
std::string fname; //!< Name of the field
};

17
src/ASM/Fields.h
View File

@ -18,7 +18,7 @@
#include <string>
class ASMbase;
class FiniteElement;
class ItgPoint;
class Vec4;
@ -73,19 +73,19 @@ public:
virtual bool valueNode(size_t node, Vector& vals) const;
//! \brief Computes the value for a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[in] x Local coordinate of evaluation point
//! \param[out] vals Values at local point in given element
virtual bool valueFE(const FiniteElement& fe, Vector& vals) const = 0;
virtual bool valueFE(const ItgPoint& x, Vector& vals) const = 0;
//! \brief Computes the value for a given global coordinate.
//! \param[in] x Global coordinate of evaluation point
//! \param[in] x Global/physical coordinate for point
//! \param[out] vals Values at given global coordinate
virtual bool valueCoor(const Vec4& x, Vector& vals) const { return false; }
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[in] x Local coordinate of evaluation point
//! \param[out] grad Gradient at local point in given element
virtual bool gradFE(const FiniteElement& fe, Matrix& grad) const = 0;
virtual bool gradFE(const ItgPoint& x, Matrix& grad) const = 0;
//! \brief Computes the gradient for a given global coordinate.
//! \param[in] x Global coordinate of evaluation point
@ -93,10 +93,9 @@ public:
virtual bool gradCoor(const Vec4& x, Matrix& grad) const { return false; }
//! \brief Computes the hessian for a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[in] x Local coordinate of evaluation point
//! \param[out] H Hessian at local point in given element
virtual bool hessianFE(const FiniteElement& fe, Matrix3D& H) const
{ return false; }
virtual bool hessianFE(const ItgPoint& x, Matrix3D& H) const { return false; }
protected:
unsigned char nf; //!< Number of field components

15
src/ASM/FiniteElement.h
View File

@ -14,6 +14,7 @@
#ifndef _FINITE_ELEMENT_H
#define _FINITE_ELEMENT_H
#include "ItgPoint.h"
#include "MatVec.h"
#include "Vec3.h"
#include "Tensor.h"
@ -23,12 +24,12 @@
\brief Class representing a finite element.
*/
class FiniteElement
class FiniteElement : public ItgPoint
{
public:
//! \brief Default constructor.
explicit FiniteElement(size_t n = 0, size_t i = 0) : iGP(i), N(n), Te(3)
{ iel = p = q = r = 0; u = v = w = xi = eta = zeta = h = 0.0; detJxW = 1.0; }
explicit FiniteElement(size_t n = 0, size_t i = 0) : ItgPoint(i), N(n), Te(3)
{ p = q = r = 0; h = 0.0; detJxW = 1.0; }
//! \brief Empty destructor.
virtual ~FiniteElement() {}
@ -68,13 +69,6 @@ protected:
public:
// Gauss point quantities
size_t iGP; //!< Global integration point counter
double u; //!< First parameter of current point
double v; //!< Second parameter of current point
double w; //!< Third parameter of current point
double xi; //!< First local coordinate of current integration point
double eta; //!< Second local coordinate of current integration point
double zeta; //!< Third local coordinate of current integration point
double detJxW; //!< Weighted determinant of the coordinate mapping
Vector N; //!< Basis function values
Matrix dNdX; //!< First derivatives (gradient) of the basis functions
@ -84,7 +78,6 @@ public:
Matrix H; //!< Hessian
// Element quantities
int iel; //!< Element identifier
short int p; //!< Polynomial order of the basis in u-direction
short int q; //!< Polynomial order of the basis in v-direction
short int r; //!< Polynomial order of the basis in r-direction

74
src/ASM/ItgPoint.h Normal file
View File

@ -0,0 +1,74 @@
// $Id$
//==============================================================================
//!
//! \file ItgPoint.h
//!
//! \date Sep 30 2019
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Integration point representation.
//!
//==============================================================================
#ifndef _ITG_POINT_H
#define _ITG_POINT_H
#include <cstddef>
/*!
\brief Class representing an integration point.
*/
class ItgPoint
{
public:
//! \brief Default constructor.
explicit ItgPoint(size_t i = 0)
{
iGP = i;
iel = -1;
u = v = w = xi = eta = zeta = 0.0;
}
//! \brief Constructor initializing the spline domain parameters.
explicit ItgPoint(double a, double b = 0.0, double c = 0.0, size_t i = 0)
{
u = a;
v = b;
w = c;
iGP = i;
iel = -1;
xi = eta = zeta = 0.0;
}
//! \brief Alternative constructor initializing the spline domain parameters.
explicit ItgPoint(const double* par, size_t i = 0)
{
u = par[0];
v = par[1];
w = par[2];
iGP = i;
iel = -1;
xi = eta = zeta = 0.0;
}
//! \brief Empty destructor.
virtual ~ItgPoint() {}
size_t iGP; //!< Global integration point counter
double u; //!< First spline parameter of the point
double v; //!< Second spline parameter of the point
double w; //!< Third spline parameter of the point
int iel; //!< Identifier of the element containing this point
double xi; //!< First local coordinate within current element
double eta; //!< Second local coordinate within current element
double zeta; //!< Third local coordinate within current element
};
#endif

52
src/ASM/LR/LRSplineField2D.C
View File

@ -11,15 +11,13 @@
//!
//==============================================================================
#include "LRSplineField2D.h"
#include "ASMu2D.h"
#include "FiniteElement.h"
#include "CoordinateMapping.h"
#include "Vec3.h"
#include "LRSpline/LRSplineSurface.h"
#include <cassert>
#include "LRSplineField2D.h"
#include "ASMu2D.h"
#include "ItgPoint.h"
#include "CoordinateMapping.h"
#include "Vec3.h"
LRSplineField2D::LRSplineField2D (const ASMu2D* patch,
@ -55,15 +53,15 @@ double LRSplineField2D::valueNode (size_t node) const
}
double LRSplineField2D::valueFE (const FiniteElement& fe) const
double LRSplineField2D::valueFE (const ItgPoint& x) const
{
if (!basis) return 0.0;
// Evaluate the basis functions at the given point
int iel = basis->getElementContaining(fe.u,fe.v);
int iel = basis->getElementContaining(x.u,x.v);
auto elm = basis->getElement(iel);
Go::BasisPtsSf spline;
basis->computeBasis(fe.u,fe.v,spline,iel);
basis->computeBasis(x.u,x.v,spline,iel);
Vector Vnod;
Vnod.reserve(elm->nBasisFunctions());
@ -77,29 +75,21 @@ double LRSplineField2D::valueFE (const FiniteElement& fe) const
double LRSplineField2D::valueCoor (const Vec4& x) const
{
if (x.u) {
FiniteElement fe;
fe.u = x.u[0];
fe.v = x.u[1];
return this->valueFE(fe);
}
assert(0);
return 0.0;
// Just produce a segmentation fault, if invoked without the parameters
return this->valueFE(ItgPoint(x.u[0],x.u[1]));
}
bool LRSplineField2D::gradFE (const FiniteElement& fe, Vector& grad) const
bool LRSplineField2D::gradFE (const ItgPoint& x, Vector& grad) const
{
if (!basis) return false;
if (!surf) return false;
// Evaluate the basis functions at the given point
int iel = surf->getElementContaining(fe.u,fe.v);
int iel = surf->getElementContaining(x.u,x.v);
auto elm = surf->getElement(iel);
Go::BasisDerivsSf spline;
basis->computeBasis(fe.u,fe.v,spline,iel);
basis->computeBasis(x.u,x.v,spline,iel);
const size_t nen = elm->nBasisFunctions();
@ -130,9 +120,9 @@ bool LRSplineField2D::gradFE (const FiniteElement& fe, Vector& grad) const
if (basis != surf)
{
// Mixed formulation, the solution uses a different basis than the geometry
int iel = basis->getElementContaining(fe.u,fe.v);
int iel = basis->getElementContaining(x.u,x.v);
auto belm = basis->getElement(iel);
basis->computeBasis(fe.u,fe.v,spline,iel);
basis->computeBasis(x.u,x.v,spline,iel);
const size_t nbf = belm->nBasisFunctions();
dNdu.resize(nbf,2);
@ -153,12 +143,12 @@ bool LRSplineField2D::gradFE (const FiniteElement& fe, Vector& grad) const
}
bool LRSplineField2D::hessianFE (const FiniteElement& fe, Matrix& H) const
bool LRSplineField2D::hessianFE (const ItgPoint& x, Matrix& H) const
{
if (!basis) return false;
if (!surf) return false;
int iel = surf->getElementContaining(fe.u,fe.v);
int iel = surf->getElementContaining(x.u,x.v);
auto elm = surf->getElement(iel);
const size_t nen = elm->nBasisFunctions();
@ -176,7 +166,7 @@ bool LRSplineField2D::hessianFE (const FiniteElement& fe, Matrix& H) const
Xnod(i, j) = (*it)->cp(j-1);
if (surf == basis) {
surf->computeBasis(fe.u,fe.v,spline2,iel);
surf->computeBasis(x.u,x.v,spline2,iel);
d2Ndu2.resize(nen,2,2);
for (size_t n = 1; n <= nen; n++) {
dNdu(n,1) = spline2.basisDerivs_u[n-1];
@ -192,7 +182,7 @@ bool LRSplineField2D::hessianFE (const FiniteElement& fe, Matrix& H) const
Vnod.push_back(values((*it)->getId()+1));
}
else {
surf->computeBasis(fe.u,fe.v,spline,iel);
surf->computeBasis(x.u,x.v,spline,iel);
for (size_t n = 1; n <= nen; n++) {
dNdu(n,1) = spline.basisDerivs_u[n-1];
dNdu(n,2) = spline.basisDerivs_v[n-1];
@ -205,9 +195,9 @@ bool LRSplineField2D::hessianFE (const FiniteElement& fe, Matrix& H) const
// 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
int iel = basis->getElementContaining(fe.u,fe.v);
int iel = basis->getElementContaining(x.u,x.v);
auto belm = basis->getElement(iel);
basis->computeBasis(fe.u,fe.v,spline2,iel);
basis->computeBasis(x.u,x.v,spline2,iel);
const size_t nbf = belm->nBasisFunctions();
dNdu.resize(nbf,2);

21
src/ASM/LR/LRSplineField2D.h
View File

@ -26,8 +26,8 @@ namespace LR {
/*!
\brief Class for LR spline-based finite element scalar fields in 2D.
\details This class implements the methods required to evaluate a 2D
lr spline scalar field at a given point in parametrical or physical coordinates.
\details This class implements the methods required to evaluate a 2D LR spline
scalar field at a given point in parametrical or physical coordinates.
*/
class LRSplineField2D : public FieldBase
@ -37,10 +37,11 @@ public:
//! \param[in] patch The spline patch on which the field is to be defined
//! \param[in] v Array of control point field values
//! \param[in] basis Basis to use from patch
//! \param[in] cmp Component to use from source field. Pass 0 to use vector as-is.
//! \param[in] cmp Component to use from source field.
//! Pass 0 to use vector as-is.
//! \param[in] name Name of spline field
LRSplineField2D(const ASMu2D* patch, const RealArray& v,
char basis = 1, char cmp = 1, const char* name = nullptr);
char basis = 1, char cmp = 1, const char* name = nullptr);
//! \brief Empty destructor.
virtual ~LRSplineField2D() {}
@ -52,22 +53,22 @@ public:
virtual double valueNode(size_t node) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
virtual double valueFE(const FiniteElement& fe) const;
//! \param[in] x Local coordinate of evaluation point
virtual double valueFE(const ItgPoint& x) const;
//! \brief Computes the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point
virtual double valueCoor(const Vec4& x) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[in] x Local coordinate of evaluation point
//! \param[out] grad Gradient of solution in a given local coordinate
virtual bool gradFE(const FiniteElement& fe, Vector& grad) const;
virtual bool gradFE(const ItgPoint& x, Vector& grad) const;
//! \brief Computes the hessian for a given local coordinate.
//! \param[in] fe Finite element quantities
//! \param[in] x Local coordinate of evaluation point
//! \param[out] H Hessian of solution in a given local coordinate
virtual bool hessianFE(const FiniteElement& fe, Matrix& H) const;
virtual bool hessianFE(const ItgPoint& x, Matrix& H) const;
protected:
const LR::LRSplineSurface* basis; //!< Spline basis description

47
src/ASM/LR/LRSplineField3D.C
View File

@ -11,16 +11,14 @@
//!
//==============================================================================
#include "LRSpline/LRSplineVolume.h"
#include "LRSplineField3D.h"
#include "ASMu3D.h"
#include "FiniteElement.h"
#include "CoordinateMapping.h"
#include "Vec3.h"
#include "LRSpline/LRSplineVolume.h"
#include <cassert>
LRSplineField3D::LRSplineField3D (const ASMu3D* patch,
const RealArray& v, char nbasis,
@ -55,15 +53,15 @@ double LRSplineField3D::valueNode (size_t node) const
}
double LRSplineField3D::valueFE (const FiniteElement& fe) const
double LRSplineField3D::valueFE (const ItgPoint& x) const
{
if (!basis) return 0.0;
// Evaluate the basis functions at the given point
int iel = basis->getElementContaining(fe.u,fe.v,fe.w);
int iel = basis->getElementContaining(x.u,x.v,x.w);
auto elm = basis->getElement(iel);
Go::BasisPts spline;
basis->computeBasis(fe.u,fe.v,fe.w,spline,iel);
basis->computeBasis(x.u,x.v,x.w,spline,iel);
Vector Vnod;
Vnod.reserve(elm->nBasisFunctions());
@ -77,30 +75,21 @@ double LRSplineField3D::valueFE (const FiniteElement& fe) const
double LRSplineField3D::valueCoor (const Vec4& x) const
{
if (x.u) {
FiniteElement fe;
fe.u = x.u[0];
fe.v = x.u[1];
fe.w = x.u[2];
return this->valueFE(fe);
}
assert(0);
return 0.0;
// Just produce a segmentation fault, if invoked without the parameters
return this->valueFE(ItgPoint(x.u[0],x.u[1],x.u[2]));
}
bool LRSplineField3D::gradFE (const FiniteElement& fe, Vector& grad) const
bool LRSplineField3D::gradFE (const ItgPoint& x, Vector& grad) const
{
if (!basis) return false;
if (!vol) return false;
// Evaluate the basis functions at the given point
int iel = vol->getElementContaining(fe.u,fe.v,fe.w);
int iel = vol->getElementContaining(x.u,x.v,x.w);
auto elm = vol->getElement(iel);
Go::BasisDerivs spline;
vol->computeBasis(fe.u,fe.v,fe.w,spline,iel);
vol->computeBasis(x.u,x.v,x.w,spline,iel);
const size_t nen = elm->nBasisFunctions();
@ -132,9 +121,9 @@ bool LRSplineField3D::gradFE (const FiniteElement& fe, Vector& grad) const
if (basis != vol)
{
// Mixed formulation, the solution uses a different basis than the geometry
int iel = basis->getElementContaining(fe.u,fe.v,fe.w);
int iel = basis->getElementContaining(x.u,x.v,x.w);
auto belm = basis->getElement(iel);
basis->computeBasis(fe.u,fe.v,fe.w,spline,iel);
basis->computeBasis(x.u,x.v,x.w,spline,iel);
const size_t nbf = belm->nBasisFunctions();
dNdu.resize(nbf,3);
@ -156,12 +145,12 @@ bool LRSplineField3D::gradFE (const FiniteElement& fe, Vector& grad) const
}
bool LRSplineField3D::hessianFE (const FiniteElement& fe, Matrix& H) const
bool LRSplineField3D::hessianFE (const ItgPoint& x, Matrix& H) const
{
if (!basis) return false;
if (!vol) return false;
int iel = vol->getElementContaining(fe.u,fe.v,fe.w);
int iel = vol->getElementContaining(x.u,x.v,x.w);
auto elm = vol->getElement(iel);
const size_t nen = elm->nBasisFunctions();
@ -179,7 +168,7 @@ bool LRSplineField3D::hessianFE (const FiniteElement& fe, Matrix& H) const
Xnod(i, j) = (*it)->cp(j-1);
if (vol == basis) {
vol->computeBasis(fe.u,fe.v,fe.w,spline2,iel);
vol->computeBasis(x.u,x.v,x.w,spline2,iel);
d2Ndu2.resize(nen,3,3);
for (size_t n = 1; n <= nen; n++) {
dNdu(n,1) = spline2.basisDerivs_u[n-1];
@ -199,7 +188,7 @@ bool LRSplineField3D::hessianFE (const FiniteElement& fe, Matrix& H) const
Vnod.push_back(values((*it)->getId()+1));
}
else {
vol->computeBasis(fe.u,fe.v,fe.w,spline,iel);
vol->computeBasis(x.u,x.v,x.w,spline,iel);
for (size_t n = 1; n <= nen; n++) {
dNdu(n,1) = spline.basisDerivs_u[n-1];
dNdu(n,2) = spline.basisDerivs_v[n-1];
@ -213,9 +202,9 @@ bool LRSplineField3D::hessianFE (const FiniteElement& fe, Matrix& H) const
// 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
int iel = basis->getElementContaining(fe.u,fe.v,fe.w);
int iel = basis->getElementContaining(x.u,x.v,x.w);
auto belm = basis->getElement(iel);
basis->computeBasis(fe.u,fe.v,fe.w,spline2,iel);
basis->computeBasis(x.u,x.v,x.w,spline2,iel);
const size_t nbf = belm->nBasisFunctions();
dNdu.resize(nbf,3);

19
src/ASM/LR/LRSplineField3D.h
View File

@ -26,8 +26,8 @@ namespace LR {
/*!
\brief Class for LR spline-based finite element scalar fields in 3D.
\details This class implements the functions required to evaluate a 3D
LR spline scalar field at a given point in parametrical or physical coordinates.
\details This class implements the methods required to evaluate a 3D LR spline
scalar field at a given point in parametrical or physical coordinates.
*/
class LRSplineField3D : public FieldBase
@ -37,7 +37,8 @@ public:
//! \param[in] patch The spline patch on which the field is to be defined
//! \param[in] v Array of control point field values
//! \param[in] basis Basis to use from patch
//! \param[in] cmp Component to use from source field. Pass 0 to use vector as-is.
//! \param[in] cmp Component to use from source field.
//! Pass 0 to use vector as-is.
//! \param[in] name Name of spline field
LRSplineField3D(const ASMu3D* patch, const RealArray& v,
char basis = 1, char cmp = 1, const char* name = nullptr);
@ -52,22 +53,22 @@ public:
virtual double valueNode(size_t node) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
virtual double valueFE(const FiniteElement& fe) const;
//! \param[in] x Local coordinate of evaluation point
virtual double valueFE(const ItgPoint& x) const;
//! \brief Computes the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point
virtual double valueCoor(const Vec4& x) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[in] x Local coordinate of evaluation point
//! \param[out] grad Gradient of solution in a given local coordinate
virtual bool gradFE(const FiniteElement& fe, Vector& grad) const;
virtual bool gradFE(const ItgPoint& x, Vector& grad) const;
//! \brief Computes the hessian for a given local coordinate.
//! \param[in] fe Finite element quantities
//! \param[in] x Local coordinate of evaluation point
//! \param[out] H Hessian of solution in a given local coordinate
virtual bool hessianFE(const FiniteElement& fe, Matrix& H) const;
virtual bool hessianFE(const ItgPoint& x, Matrix& H) const;
protected:
const LR::LRSplineVolume* basis; //!< Spline basis description

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

@ -11,15 +11,13 @@
//!
//==============================================================================
#include "LRSplineFields2D.h"
#include "ASMu2D.h"
#include "FiniteElement.h"
#include "CoordinateMapping.h"
#include "Vec3.h"
#include "LRSpline/LRSplineSurface.h"
#include <cassert>
#include "LRSplineFields2D.h"
#include "ASMu2D.h"
#include "ItgPoint.h"
#include "CoordinateMapping.h"
#include "Vec3.h"
LRSplineFields2D::LRSplineFields2D (const ASMu2D* patch,
@ -72,16 +70,16 @@ bool LRSplineFields2D::valueNode (size_t node, Vector& vals) const
}
bool LRSplineFields2D::valueFE (const FiniteElement& fe, Vector& vals) const
bool LRSplineFields2D::valueFE (const ItgPoint& x, Vector& vals) const
{
if (!basis) return false;
// Evaluate the basis functions at the given point
int iel = basis->getElementContaining(fe.u,fe.v);
int iel = basis->getElementContaining(x.u,x.v);
auto elm = basis->getElement(iel);
Go::BasisPtsSf spline;
basis->computeBasis(fe.u,fe.v,spline,iel);
basis->computeBasis(x.u,x.v,spline,iel);
// Evaluate the solution field at the given point
Matrix Vnod(nf, elm->nBasisFunctions());
@ -99,29 +97,23 @@ bool LRSplineFields2D::valueFE (const FiniteElement& fe, Vector& vals) const
bool LRSplineFields2D::valueCoor (const Vec4& x, Vector& vals) const
{
if (x.u) {
FiniteElement fe;
fe.u = x.u[0];
fe.v = x.u[1];
if (x.u)
return this->valueFE(ItgPoint(x.u[0],x.u[1]),vals);
return this->valueFE(fe, vals);
}
assert(0);
return false;
}
bool LRSplineFields2D::gradFE (const FiniteElement& fe, Matrix& grad) const
bool LRSplineFields2D::gradFE (const ItgPoint& x, Matrix& grad) const
{
if (!basis) return false;
if (!surf) return false;
// Evaluate the basis functions at the given point
int iel = surf->getElementContaining(fe.u,fe.v);
int iel = surf->getElementContaining(x.u,x.v);
auto elm = surf->getElement(iel);
Go::BasisDerivsSf spline;
surf->computeBasis(fe.u,fe.v,spline,iel);
surf->computeBasis(x.u,x.v,spline,iel);
const size_t nen = elm->nBasisFunctions();
@ -147,9 +139,9 @@ bool LRSplineFields2D::gradFE (const FiniteElement& fe, Matrix& grad) const
if (basis != surf)
{
// Mixed formulation, the solution uses a different basis than the geometry
int iel = basis->getElementContaining(fe.u,fe.v);
int iel = basis->getElementContaining(x.u,x.v);
auto belm = basis->getElement(iel);
basis->computeBasis(fe.u,fe.v,spline,iel);
basis->computeBasis(x.u,x.v,spline,iel);
const size_t nbf = belm->nBasisFunctions();
dNdu.resize(nbf,2);
@ -182,12 +174,12 @@ bool LRSplineFields2D::gradFE (const FiniteElement& fe, Matrix& grad) const
}
bool LRSplineFields2D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
bool LRSplineFields2D::hessianFE (const ItgPoint& x, Matrix3D& H) const
{
if (!basis) return false;
if (!surf) return false;
int iel = surf->getElementContaining(fe.u,fe.v);
int iel = surf->getElementContaining(x.u,x.v);
auto elm = surf->getElement(iel);
const size_t nen = elm->nBasisFunctions();
@ -204,7 +196,7 @@ bool LRSplineFields2D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
Xnod(j,i) = (*it)->cp(j-1);
if (surf == basis) {
surf->computeBasis(fe.u,fe.v,spline2,iel);
surf->computeBasis(x.u,x.v,spline2,iel);
dNdu.resize(nen,2);
d2Ndu2.resize(nen,2,2);
@ -224,7 +216,7 @@ bool LRSplineFields2D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
Vnod(j,i) = values((*it)->getId()*nf+j);
}
else {
surf->computeBasis(fe.u,fe.v,spline,iel);
surf->computeBasis(x.u,x.v,spline,iel);
dNdu.resize(nen,2);
for (size_t n = 1; n <= nen; n++) {
@ -239,9 +231,9 @@ bool LRSplineFields2D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
// 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
int iel = basis->getElementContaining(fe.u,fe.v);
int iel = basis->getElementContaining(x.u,x.v);
auto belm = basis->getElement(iel);
basis->computeBasis(fe.u,fe.v,spline2,iel);
basis->computeBasis(x.u,x.v,spline2,iel);
const size_t nbf = belm->nBasisFunctions();
dNdu.resize(nbf,2);

16
src/ASM/LR/LRSplineFields2D.h
View File

@ -26,8 +26,8 @@ namespace LR {
/*!
\brief Class for LR spline-based finite element vector fields in 2D.
\details This class implements the methods required to evaluate a 2D
LR spline vector field at a given point in parametrical or physical coordinates.
\details This class implements the methods required to evaluate a 2D LR spline
vector field at a given point in parametrical or physical coordinates.
*/
class LRSplineFields2D : public Fields
@ -62,9 +62,9 @@ public:
bool valueNode(size_t node, Vector& vals) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[in] x Local coordinate of evaluation point
//! \param[out] vals Values in local point in given element
bool valueFE(const FiniteElement& fe, Vector& vals) const;
bool valueFE(const ItgPoint& x, Vector& vals) const;
//! \brief Computes the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point
@ -72,14 +72,14 @@ public:
bool valueCoor(const Vec4& x, Vector& vals) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[in] x Local coordinate of evaluation point
//! \param[out] grad Gradient of solution in a given local coordinate
bool gradFE(const FiniteElement& fe, Matrix& grad) const;
bool gradFE(const ItgPoint& x, Matrix& grad) const;
//! \brief Computes the hessian for a given local coordinate.
//! \param[in] fe Finite element quantities
//! \param[in] x Local coordinate of evaluation point
//! \param[out] H Hessian of solution in a given local coordinate
virtual bool hessianFE(const FiniteElement& fe, Matrix3D& H) const;
virtual bool hessianFE(const ItgPoint& x, Matrix3D& H) const;
protected:
const LR::LRSplineSurface* basis; //!< Spline basis description

22
src/ASM/LR/LRSplineFields2Dmx.C
View File

@ -11,14 +11,14 @@
//!
//==============================================================================
#include "LRSpline/LRSplineSurface.h"
#include "LRSplineFields2Dmx.h"
#include "ASMu2Dmx.h"
#include "FiniteElement.h"
#include "ItgPoint.h"
#include "CoordinateMapping.h"
#include "Utilities.h"
#include "LRSpline/LRSplineSurface.h"
LRSplineFields2Dmx::LRSplineFields2Dmx (const ASMu2Dmx* patch,
const RealArray& v, char basis,
@ -48,7 +48,7 @@ bool LRSplineFields2Dmx::valueNode (size_t node, Vector& vals) const
}
bool LRSplineFields2Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
bool LRSplineFields2Dmx::valueFE (const ItgPoint& x, Vector& vals) const
{
if (!surf) return false;
@ -60,10 +60,10 @@ bool LRSplineFields2Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
for (const auto& it : bases) {
const LR::LRSplineSurface* basis = surf->getBasis(it);
int iel = basis->getElementContaining(fe.u,fe.v);
int iel = basis->getElementContaining(x.u,x.v);
auto elm = basis->getElement(iel);
Go::BasisPtsSf spline;
basis->computeBasis(fe.u,fe.v,spline,iel);
basis->computeBasis(x.u,x.v,spline,iel);
// Evaluate the solution field at the given point
@ -81,16 +81,16 @@ bool LRSplineFields2Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
}
bool LRSplineFields2Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
bool LRSplineFields2Dmx::gradFE (const ItgPoint& x, Matrix& grad) const
{
if (!surf) return false;
// Evaluate the basis functions at the given point
Go::BasisDerivsSf spline;
const LR::LRSplineSurface* gsurf = surf->getBasis(ASMmxBase::geoBasis);
int iel = gsurf->getElementContaining(fe.u,fe.v);
int iel = gsurf->getElementContaining(x.u,x.v);
auto elm = gsurf->getElement(iel);
gsurf->computeBasis(fe.u,fe.v,spline,iel);
gsurf->computeBasis(x.u,x.v,spline,iel);
const size_t nen = elm->nBasisFunctions();
@ -116,9 +116,9 @@ bool LRSplineFields2Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
size_t row = 1;
for (const auto& it : bases) {
const LR::LRSplineSurface* basis = surf->getBasis(it);
int iel = basis->getElementContaining(fe.u,fe.v);
int iel = basis->getElementContaining(x.u,x.v);
auto belm = basis->getElement(iel);
basis->computeBasis(fe.u,fe.v,spline,iel);
basis->computeBasis(x.u,x.v,spline,iel);
const size_t nbf = belm->nBasisFunctions();
dNdu.resize(nbf,2);

14
src/ASM/LR/LRSplineFields2Dmx.h
View File

@ -27,8 +27,8 @@ namespace LR {
/*!
\brief Class for mixed LR spline-based finite element vector fields in 2D.
\details This class implements the methods required to evaluate a 2D
mixed LR spline vector field at a given point in parametrical or physical coordinates.
\details This class implements the methods required to evaluate a 2D LR spline
mixed vector field at a given point in parametrical or physical coordinates.
*/
class LRSplineFields2Dmx : public Fields
@ -53,18 +53,18 @@ public:
bool valueNode(size_t node, Vector& vals) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[in] x Local coordinate of evaluation point
//! \param[out] vals Values in local point in given element
bool valueFE(const FiniteElement& fe, Vector& vals) const;
bool valueFE(const ItgPoint& x, Vector& vals) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[in] x Local coordinate of evaluation point
//! \param[out] grad Gradient of solution in a given local coordinate
bool gradFE(const FiniteElement& fe, Matrix& grad) const;
bool gradFE(const ItgPoint& x, Matrix& grad) const;
protected:
const ASMu2Dmx* surf; //!< Patch description
std::set<int> bases; //!< Bases to use
std::set<int> bases; //!< Bases to use
};
#endif

51
src/ASM/LR/LRSplineFields3D.C
View File

@ -11,15 +11,13 @@
//!
//==============================================================================
#include "LRSplineFields3D.h"
#include "ASMu3D.h"
#include "FiniteElement.h"
#include "CoordinateMapping.h"
#include "Vec3.h"
#include "LRSpline/LRSplineVolume.h"
#include <cassert>
#include "LRSplineFields3D.h"
#include "ASMu3D.h"
#include "ItgPoint.h"
#include "CoordinateMapping.h"
#include "Vec3.h"
LRSplineFields3D::LRSplineFields3D (const ASMu3D* patch,
@ -72,16 +70,16 @@ bool LRSplineFields3D::valueNode (size_t node, Vector& vals) const
}
bool LRSplineFields3D::valueFE (const FiniteElement& fe, Vector& vals) const
bool LRSplineFields3D::valueFE (const ItgPoint& x, Vector& vals) const
{
if (!basis) return false;
// Evaluate the basis functions at the given point
int iel = basis->getElementContaining(fe.u,fe.v,fe.w);
int iel = basis->getElementContaining(x.u,x.v,x.w);
auto elm = basis->getElement(iel);
Go::BasisPts spline;
basis->computeBasis(fe.u,fe.v,fe.w,spline,iel);
basis->computeBasis(x.u,x.v,x.w,spline,iel);
// Evaluate the solution field at the given point
Matrix Vnod(nf, elm->nBasisFunctions());
@ -99,30 +97,23 @@ bool LRSplineFields3D::valueFE (const FiniteElement& fe, Vector& vals) const
bool LRSplineFields3D::valueCoor (const Vec4& x, Vector& vals) const
{
if (x.u) {
FiniteElement fe;
fe.u = x.u[0];
fe.v = x.u[1];
fe.w = x.u[2];
if (x.u)
return this->valueFE(ItgPoint(x.u[0],x.u[1],x.u[2]),vals);
return this->valueFE(fe, vals);
}
assert(0);
return false;
}
bool LRSplineFields3D::gradFE (const FiniteElement& fe, Matrix& grad) const
bool LRSplineFields3D::gradFE (const ItgPoint& x, Matrix& grad) const
{
if (!basis) return false;
if (!vol) return false;
// Evaluate the basis functions at the given point
int iel = vol->getElementContaining(fe.u,fe.v,fe.w);
int iel = vol->getElementContaining(x.u,x.v,x.w);
auto elm = vol->getElement(iel);
Go::BasisDerivs spline;
vol->computeBasis(fe.u,fe.v,fe.w,spline,iel);
vol->computeBasis(x.u,x.v,x.w,spline,iel);
const size_t nen = elm->nBasisFunctions();
@ -149,9 +140,9 @@ bool LRSplineFields3D::gradFE (const FiniteElement& fe, Matrix& grad) const
if (basis != vol)
{
// Mixed formulation, the solution uses a different basis than the geometry
int iel = basis->getElementContaining(fe.u,fe.v,fe.w);
int iel = basis->getElementContaining(x.u,x.v,x.w);
auto belm = basis->getElement(iel);
basis->computeBasis(fe.u,fe.v,fe.w,spline,iel);
basis->computeBasis(x.u,x.v,x.w,spline,iel);
const size_t nbf = belm->nBasisFunctions();
dNdu.resize(nbf,3);
@ -185,12 +176,12 @@ bool LRSplineFields3D::gradFE (const FiniteElement& fe, Matrix& grad) const
}
bool LRSplineFields3D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
bool LRSplineFields3D::hessianFE (const ItgPoint& x, Matrix3D& H) const
{
if (!basis) return false;
if (!vol) return false;
int iel = vol->getElementContaining(fe.u,fe.v,fe.w);
int iel = vol->getElementContaining(x.u,x.v,x.w);
auto elm = vol->getElement(iel);
const size_t nen = elm->nBasisFunctions();
@ -207,7 +198,7 @@ bool LRSplineFields3D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
Xnod(j,i) = (*it)->cp(j-1);
if (vol == basis) {
vol->computeBasis(fe.u,fe.v,fe.w,spline2,iel);
vol->computeBasis(x.u,x.v,x.w,spline2,iel);
dNdu.resize(nen,3);
d2Ndu2.resize(nen,3,3);
@ -231,7 +222,7 @@ bool LRSplineFields3D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
Vnod(j,i) = values((*it)->getId()*nf+j);
}
else {
vol->computeBasis(fe.u,fe.v,fe.w,spline,iel);
vol->computeBasis(x.u,x.v,x.w,spline,iel);
dNdu.resize(nen,3);
for (size_t n = 1; n <= nen; n++) {
@ -247,9 +238,9 @@ bool LRSplineFields3D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
// 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
int iel = basis->getElementContaining(fe.u,fe.v,fe.w);
int iel = basis->getElementContaining(x.u,x.v,x.w);
auto belm = basis->getElement(iel);
basis->computeBasis(fe.u,fe.v,fe.w,spline2,iel);
basis->computeBasis(x.u,x.v,x.w,spline2,iel);
const size_t nbf = belm->nBasisFunctions();
dNdu.resize(nbf,3);

16
src/ASM/LR/LRSplineFields3D.h
View File

@ -26,8 +26,8 @@ namespace LR {
/*!
\brief Class for LR spline-based finite element vector fields in 3D.
\details This class implements the methods required to evaluate a 3D
LR spline vector field at a given point in parametrical or physical coordinates.
\details This class implements the methods required to evaluate a 3D LR spline
vector field at a given point in parametrical or physical coordinates.
*/
class LRSplineFields3D : public Fields
@ -62,9 +62,9 @@ public:
bool valueNode(size_t node, Vector& vals) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[in] x Local coordinate of evaluation point
//! \param[out] vals Values in local point in given element
bool valueFE(const FiniteElement& fe, Vector& vals) const;
bool valueFE(const ItgPoint& x, Vector& vals) const;
//! \brief Computes the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point
@ -72,14 +72,14 @@ public:
bool valueCoor(const Vec4& x, Vector& vals) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[in] x Local coordinate of evaluation point
//! \param[out] grad Gradient of solution in a given local coordinate
bool gradFE(const FiniteElement& fe, Matrix& grad) const;
bool gradFE(const ItgPoint& x, Matrix& grad) const;
//! \brief Computes the hessian for a given local coordinate.
//! \param[in] fe Finite element quantities
//! \param[in] x Local coordinate of evaluation point
//! \param[out] H Hessian of solution in a given local coordinate
virtual bool hessianFE(const FiniteElement& fe, Matrix3D& H) const;
virtual bool hessianFE(const ItgPoint& x, Matrix3D& H) const;
protected:
const LR::LRSplineVolume* basis; //!< Spline basis description

22
src/ASM/LR/LRSplineFields3Dmx.C
View File

@ -11,14 +11,14 @@
//!
//==============================================================================
#include "LRSpline/LRSplineVolume.h"
#include "LRSplineFields3Dmx.h"
#include "ASMu3Dmx.h"
#include "FiniteElement.h"
#include "ItgPoint.h"
#include "CoordinateMapping.h"
#include "Utilities.h"
#include "LRSpline/LRSplineVolume.h"
LRSplineFields3Dmx::LRSplineFields3Dmx (const ASMu3Dmx* patch,
const RealArray& v, char basis,
@ -48,7 +48,7 @@ bool LRSplineFields3Dmx::valueNode (size_t node, Vector& vals) const
}
bool LRSplineFields3Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
bool LRSplineFields3Dmx::valueFE (const ItgPoint& x, Vector& vals) const
{
if (!vol) return false;
@ -60,10 +60,10 @@ bool LRSplineFields3Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
for (const auto& it : bases) {
const LR::LRSplineVolume* basis = vol->getBasis(it);
int iel = basis->getElementContaining(fe.u,fe.v,fe.w);
int iel = basis->getElementContaining(x.u,x.v,x.w);
auto elm = basis->getElement(iel);
Go::BasisPts spline;
basis->computeBasis(fe.u,fe.v,fe.w,spline,iel);
basis->computeBasis(x.u,x.v,x.w,spline,iel);
// Evaluate the solution field at the given point
@ -81,16 +81,16 @@ bool LRSplineFields3Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
}
bool LRSplineFields3Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
bool LRSplineFields3Dmx::gradFE (const ItgPoint& x, Matrix& grad) const
{
if (!vol) return false;
// Evaluate the basis functions at the given point
Go::BasisDerivs spline;
const LR::LRSplineVolume* gvol = vol->getBasis(ASMmxBase::geoBasis);
int iel = gvol->getElementContaining(fe.u,fe.v,fe.w);
int iel = gvol->getElementContaining(x.u,x.v,x.w);
auto elm = gvol->getElement(iel);
gvol->computeBasis(fe.u,fe.v,fe.w,spline,iel);
gvol->computeBasis(x.u,x.v,x.w,spline,iel);
const size_t nen = elm->nBasisFunctions();
Matrix dNdu(nen,3), dNdX;
@ -116,9 +116,9 @@ bool LRSplineFields3Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
size_t row = 1;
for (const auto& it : bases) {
const LR::LRSplineVolume* basis = vol->getBasis(it);
int iel = basis->getElementContaining(fe.u,fe.v,fe.w);
int iel = basis->getElementContaining(x.u,x.v,x.w);
auto belm = basis->getElement(iel);
basis->computeBasis(fe.u,fe.v,fe.w,spline,iel);
basis->computeBasis(x.u,x.v,x.w,spline,iel);
const size_t nbf = belm->nBasisFunctions();
dNdu.resize(nbf,3);

12
src/ASM/LR/LRSplineFields3Dmx.h
View File

@ -27,8 +27,8 @@ namespace LR {
/*!
\brief Class for LR spline-based finite element vector fields in 3D.
\details This class implements the methods required to evaluate a 3D
LR spline vector field at a given point in parametrical or physical coordinates.
\details This class implements the methods required to evaluate a 3D LR spline
mixed vector field at a given point in parametrical or physical coordinates.
*/
class LRSplineFields3Dmx : public Fields
@ -53,14 +53,14 @@ public:
bool valueNode(size_t node, Vector& vals) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[in] x Local coordinate of evaluation point
//! \param[out] vals Values in local point in given element
bool valueFE(const FiniteElement& fe, Vector& vals) const;
bool valueFE(const ItgPoint& x, Vector& vals) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[in] x Local coordinate of evaluation point
//! \param[out] grad Gradient of solution in a given local coordinate
bool gradFE(const FiniteElement& fe, Matrix& grad) const;
bool gradFE(const ItgPoint& x, Matrix& grad) const;
protected:
const ASMu3Dmx* vol; //!< Patch description

40
src/ASM/LagrangeField2D.C
View File

@ -13,15 +13,14 @@
#include "LagrangeField2D.h"
#include "ASMs2DLag.h"
#include "FiniteElement.h"
#include "ItgPoint.h"
#include "Lagrange.h"
#include "CoordinateMapping.h"
LagrangeField2D::LagrangeField2D (const ASMs2DLag* patch,
const RealArray& v,
char basis, char,
const char* name) : FieldBase(name)
const RealArray& v, char, char,
const char* name) : FieldBase(name)
{
patch->getNodalCoordinates(coord);
patch->getSize(n1,n2);
@ -42,14 +41,19 @@ double LagrangeField2D::valueNode (size_t node) const
}
double LagrangeField2D::valueFE (const FiniteElement& fe) const
double LagrangeField2D::valueFE (const ItgPoint& x) const
{
if (x.iel < 1 || (size_t)x.iel > nelm)
{
std::cerr <<" *** LagrangeField2D::valueFE: Element index "<< x.iel
<<" out of range [1,"<< nelm <<"]."<< std::endl;
return 0.0;
}
Vector N;
Lagrange::computeBasis(N,p1,fe.xi,p2,fe.eta);
Lagrange::computeBasis(N,p1,x.xi,p2,x.eta);
const int nel1 = (n1-1)/p1;
div_t divresult = div(fe.iel,nel1);
div_t divresult = div(x.iel,(n1-1)/p1);
int iel1 = divresult.rem;
int iel2 = divresult.quot;
const int node1 = p1*iel1-1;
@ -68,18 +72,22 @@ double LagrangeField2D::valueFE (const FiniteElement& fe) const
}
bool LagrangeField2D::gradFE (const FiniteElement& fe, Vector& grad) const
bool LagrangeField2D::gradFE (const ItgPoint& x, Vector& grad) const
{
grad.resize(2,true);
if (x.iel < 1 || (size_t)x.iel > nelm)
{
std::cerr <<" *** LagrangeField2D::gradFE: Element index "<< x.iel
<<" out of range [1,"<< nelm <<"]."<< std::endl;
return false;
}
Vector Vnod;
Matrix dNdu;
if (!Lagrange::computeBasis(Vnod,dNdu,p1,fe.xi,p2,fe.eta))
if (!Lagrange::computeBasis(Vnod,dNdu,p1,x.xi,p2,x.eta))
return false;
const int nel1 = (n1-1)/p1;
div_t divresult = div(fe.iel,nel1);
div_t divresult = div(x.iel,(n1-1)/p1);
int iel1 = divresult.rem;
int iel2 = divresult.quot;
const int node1 = p1*iel1-1;
@ -98,6 +106,8 @@ bool LagrangeField2D::gradFE (const FiniteElement& fe, Vector& grad) const
}
Matrix Jac, dNdX;
utl::Jacobian(Jac,dNdX,Xnod,dNdu,false);
if (!utl::Jacobian(Jac,dNdX,Xnod,dNdu))
return false; // Singular Jacobian
return dNdX.multiply(Vnod,grad);
}

8
src/ASM/LagrangeField2D.h
View File

@ -48,13 +48,13 @@ public:
double valueNode(size_t node) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
double valueFE(const FiniteElement& fe) const;
//! \param[in] x Local coordinate of evaluation point
double valueFE(const ItgPoint& x) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[in] x Local coordinate of evaluation point
//! \param[out] grad Gradient of solution in a given local coordinate
bool gradFE(const FiniteElement& fe, Vector& grad) const;
bool gradFE(const ItgPoint& x, Vector& grad) const;
protected:
Matrix coord; //!< Matrix of nodal coordinates

34
src/ASM/LagrangeField3D.C
View File

@ -13,14 +13,13 @@
#include "LagrangeField3D.h"
#include "ASMs3DLag.h"
#include "FiniteElement.h"
#include "ItgPoint.h"
#include "Lagrange.h"
#include "CoordinateMapping.h"
LagrangeField3D::LagrangeField3D (const ASMs3DLag* patch,
const RealArray& v,
char basis, char,
const RealArray& v, char, char,
const char* name) : FieldBase(name)
{
patch->getNodalCoordinates(coord);
@ -42,15 +41,22 @@ double LagrangeField3D::valueNode (size_t node) const
}
double LagrangeField3D::valueFE (const FiniteElement& fe) const
double LagrangeField3D::valueFE (const ItgPoint& x) const
{
if (x.iel < 1 || (size_t)x.iel > nelm)
{
std::cerr <<" *** LagrangeField3D::valueFE: Element index "<< x.iel
<<" out of range [1,"<< nelm <<"]."<< std::endl;
return 0.0;
}
Vector N;
Lagrange::computeBasis(N,p1,fe.xi,p2,fe.eta,p3,fe.zeta);
Lagrange::computeBasis(N,p1,x.xi,p2,x.eta,p3,x.zeta);
const int nel1 = (n1-1)/p1;
const int nel2 = (n2-1)/p2;
div_t divresult = div(fe.iel,nel1*nel2);
div_t divresult = div(x.iel,nel1*nel2);
int iel2 = divresult.rem;
int iel3 = divresult.quot;
divresult = div(iel2,nel1);
@ -74,19 +80,25 @@ double LagrangeField3D::valueFE (const FiniteElement& fe) const
}
bool LagrangeField3D::gradFE (const FiniteElement& fe, Vector& grad) const
bool LagrangeField3D::gradFE (const ItgPoint& x, Vector& grad) const
{
grad.resize(3,true);
if (x.iel < 1 || (size_t)x.iel > nelm)
{
std::cerr <<" *** LagrangeField3D::gradFE: Element index "<< x.iel
<<" out of range [1,"<< nelm <<"]."<< std::endl;
return false;
}
Vector Vnod;
Matrix dNdu;
if (!Lagrange::computeBasis(Vnod,dNdu,p1,fe.xi,p2,fe.eta,p3,fe.zeta))
if (!Lagrange::computeBasis(Vnod,dNdu,p1,x.xi,p2,x.eta,p3,x.zeta))
return false;
const int nel1 = (n1-1)/p1;
const int nel2 = (n2-1)/p2;
div_t divresult = div(fe.iel,nel1*nel2);
div_t divresult = div(x.iel,nel1*nel2);
int iel2 = divresult.rem;
int iel3 = divresult.quot;
divresult = div(iel2,nel1);
@ -110,6 +122,8 @@ bool LagrangeField3D::gradFE (const FiniteElement& fe, Vector& grad) const
}
Matrix Jac, dNdX;
utl::Jacobian(Jac,dNdX,Xnod,dNdu);
if (!utl::Jacobian(Jac,dNdX,Xnod,dNdu))
return false; // Singular Jacobian
return dNdX.multiply(Vnod,grad);
}

8
src/ASM/LagrangeField3D.h
View File

@ -48,13 +48,13 @@ public:
double valueNode(size_t node) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
double valueFE(const FiniteElement& fe) const;
//! \param[in] x Local coordinate of evaluation point
double valueFE(const ItgPoint& x) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[in] x Local coordinate of evaluation point
//! \param[out] grad Gradient of solution in a given local coordinate
bool gradFE(const FiniteElement& fe, Vector& grad) const;
bool gradFE(const ItgPoint& x, Vector& grad) const;
protected:
Matrix coord; //!< Matrix of nodel coordinates

43
src/ASM/LagrangeFields2D.C
View File

@ -13,15 +13,14 @@
#include "LagrangeFields2D.h"
#include "ASMs2DLag.h"
#include "FiniteElement.h"
#include "ItgPoint.h"
#include "Lagrange.h"
#include "CoordinateMapping.h"
LagrangeFields2D::LagrangeFields2D (const ASMs2DLag* patch,
const RealArray& v,
char basis,
const char* name) : Fields(name)
const RealArray& v, char,
const char* name) : Fields(name)
{
patch->getNodalCoordinates(coord);
patch->getSize(n1,n2);
@ -48,17 +47,21 @@ bool LagrangeFields2D::valueNode (size_t node, Vector& vals) const
}
bool LagrangeFields2D::valueFE (const FiniteElement& fe, Vector& vals) const
bool LagrangeFields2D::valueFE (const ItgPoint& x, Vector& vals) const
{
vals.resize(nf,true);
if (x.iel < 1 || (size_t)x.iel > nelm)
{
std::cerr <<" *** LagrangeFields2D::valueFE: Element index "<< x.iel
<<" out of range [1,"<< nelm <<"]."<< std::endl;
return false;
}
Vector N;
if (!Lagrange::computeBasis(N,p1,fe.xi,p2,fe.eta))
if (!Lagrange::computeBasis(N,p1,x.xi,p2,x.eta))
return false;
const int nel1 = (n1-1)/p1;
div_t divresult = div(fe.iel,nel1);
div_t divresult = div(x.iel,(n1-1)/p1);
int iel1 = divresult.rem;
int iel2 = divresult.quot;
const int node1 = p1*iel1-1;
@ -71,25 +74,29 @@ bool LagrangeFields2D::valueFE (const FiniteElement& fe, Vector& vals) const
int dof = nf*(n1*(j-1) + i-1) + 1;
double value = N(locNode);
for (int k = 1; k <= nf; k++, dof++)
vals(k) += values(dof)*value;
vals(k) += values(dof)*value;
}
return true;
}
bool LagrangeFields2D::gradFE (const FiniteElement& fe, Matrix& grad) const
bool LagrangeFields2D::gradFE (const ItgPoint& x, Matrix& grad) const
{
grad.resize(nf,2,true);
if (x.iel < 1 || (size_t)x.iel > nelm)
{
std::cerr <<" *** LagrangeFields2D::gradFE: Element index "<< x.iel
<<" out of range [1,"<< nelm <<"]."<< std::endl;
return false;
}
Vector N;
Matrix dNdu;
if (!Lagrange::computeBasis(N,dNdu,p1,fe.xi,p2,fe.eta))
if (!Lagrange::computeBasis(N,dNdu,p1,x.xi,p2,x.eta))
return false;
const int nel1 = (n1-1)/p1;
div_t divresult = div(fe.iel,nel1);
div_t divresult = div(x.iel,(n1-1)/p1);
int iel1 = divresult.rem;
int iel2 = divresult.quot;
const int node1 = p1*iel1-1;
@ -108,8 +115,8 @@ bool LagrangeFields2D::gradFE (const FiniteElement& fe, Matrix& grad) const
}
Matrix Jac, dNdX;
utl::Jacobian(Jac,dNdX,Xnod,dNdu);
grad.multiply(Vnod,dNdX);
if (!utl::Jacobian(Jac,dNdX,Xnod,dNdu))
return false; // Singular Jacobian
return true;
return !grad.multiply(Vnod,dNdX).empty();
}

8
src/ASM/LagrangeFields2D.h
View File

@ -48,14 +48,14 @@ public:
bool valueNode(size_t node, Vector& vals) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[in] x Local coordinate of evaluation point
//! \param[out] vals Values in local point in given element
bool valueFE(const FiniteElement& fe, Vector& vals) const;
bool valueFE(const ItgPoint& x, Vector& vals) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[in] x Local coordinate of evaluation point
//! \param[out] grad Gradient of solution in a given local coordinate
bool gradFE(const FiniteElement& fe, Matrix& grad) const;
bool gradFE(const ItgPoint& x, Matrix& grad) const;
protected:
Matrix coord; //!< Matrix of nodal coordinates

37
src/ASM/LagrangeFields3D.C
View File

@ -13,15 +13,14 @@
#include "LagrangeFields3D.h"
#include "ASMs3DLag.h"
#include "FiniteElement.h"
#include "ItgPoint.h"
#include "Lagrange.h"
#include "CoordinateMapping.h"
LagrangeFields3D::LagrangeFields3D (const ASMs3DLag* patch,
const RealArray& v,
char basis,
const char* name) : Fields(name)
const RealArray& v, char,
const char* name) : Fields(name)
{
patch->getNodalCoordinates(coord);
patch->getSize(n1,n2,n3);
@ -48,18 +47,24 @@ bool LagrangeFields3D::valueNode (size_t node, Vector& vals) const
}
bool LagrangeFields3D::valueFE (const FiniteElement& fe, Vector& vals) const
bool LagrangeFields3D::valueFE (const ItgPoint& x, Vector& vals) const
{
vals.resize(nf,true);
if (x.iel < 1 || (size_t)x.iel > nelm)
{
std::cerr <<" *** LagrangeFields3D::valueFE: Element index "<< x.iel
<<" out of range [1,"<< nelm <<"]."<< std::endl;
return false;
}
Vector N;
if (!Lagrange::computeBasis(N,p1,fe.xi,p2,fe.eta,p3,fe.zeta))
if (!Lagrange::computeBasis(N,p1,x.xi,p2,x.eta,p3,x.zeta))
return false;
const int nel1 = (n1-1)/p1;
const int nel2 = (n2-1)/p2;
div_t divresult = div(fe.iel,nel1*nel2);
div_t divresult = div(x.iel,nel1*nel2);
int iel2 = divresult.rem;
int iel3 = divresult.quot;
divresult = div(iel2,nel1);
@ -84,19 +89,25 @@ bool LagrangeFields3D::valueFE (const FiniteElement& fe, Vector& vals) const
}
bool LagrangeFields3D::gradFE (const FiniteElement& fe, Matrix& grad) const
bool LagrangeFields3D::gradFE (const ItgPoint& x, Matrix& grad) const
{
grad.resize(nf,3,true);
if (x.iel < 1 || (size_t)x.iel > nelm)
{
std::cerr <<" *** LagrangeFields3D::gradFE: Element index "<< x.iel
<<" out of range [1,"<< nelm <<"]."<< std::endl;
return false;
}
Vector N;
Matrix dNdu;
if (!Lagrange::computeBasis(N,dNdu,p1,fe.xi,p2,fe.eta,p3,fe.zeta))
if (!Lagrange::computeBasis(N,dNdu,p1,x.xi,p2,x.eta,p3,x.zeta))
return false;
const int nel1 = (n1-1)/p1;
const int nel2 = (n2-1)/p2;
div_t divresult = div(fe.iel,nel1*nel2);
div_t divresult = div(x.iel,nel1*nel2);
int iel2 = divresult.rem;
int iel3 = divresult.quot;
divresult = div(iel2,nel1);
@ -120,8 +131,8 @@ bool LagrangeFields3D::gradFE (const FiniteElement& fe, Matrix& grad) const
}
Matrix Jac, dNdX;
utl::Jacobian(Jac,dNdX,Xnod,dNdu);
grad.multiply(Vnod,dNdX);
if (!utl::Jacobian(Jac,dNdX,Xnod,dNdu))
return false; // Singular Jacobian
return true;
return !grad.multiply(Vnod,dNdX).empty();
}

8
src/ASM/LagrangeFields3D.h
View File

@ -48,14 +48,14 @@ public:
bool valueNode(size_t node, Vector& vals) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[in] x Local coordinate of evaluation point
//! \param[out] vals Values in local point in given element
bool valueFE(const FiniteElement& fe, Vector& vals) const;
bool valueFE(const ItgPoint& x, Vector& vals) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[in] x Local coordinate of evaluation point
//! \param[out] grad Gradient of solution in a given local coordinate
bool gradFE(const FiniteElement& fe, Matrix& grad) const;
bool gradFE(const ItgPoint& x, Matrix& grad) const;
protected:
Matrix coord; //!< Matrix of nodel coordinates

48
src/ASM/SplineField2D.C
View File

@ -15,7 +15,7 @@
#include "SplineField2D.h"
#include "ASMs2D.h"
#include "FiniteElement.h"
#include "ItgPoint.h"
#include "CoordinateMapping.h"
#include "Utilities.h"
#include "Vec3.h"
@ -60,14 +60,14 @@ double SplineField2D::valueNode (size_t node) const
}
double SplineField2D::valueFE (const FiniteElement& fe) const
double SplineField2D::valueFE (const ItgPoint& x) const
{
if (!basis) return false;
// Evaluate the basis functions at the given point
Go::BasisPtsSf spline;
#pragma omp critical
basis->computeBasis(fe.u,fe.v,spline);
basis->computeBasis(x.u,x.v,spline);
// Evaluate the solution field at the given point
IntVec ip;
@ -83,25 +83,16 @@ double SplineField2D::valueFE (const FiniteElement& fe) const
double SplineField2D::valueCoor (const Vec4& x) const
{
FiniteElement fe;
if (x.u) {
fe.u = x.u[0];
fe.v = x.u[1];
}
else {
// use with caution
Go::Point pt(3), clopt(3);
pt[0] = x[0];
pt[1] = x[1];
pt[2] = x[2];
double clo_u, clo_v, dist;
#pragma omp critical
surf->closestPoint(pt, clo_u, clo_v, clopt, dist, 1e-10);
fe.u = clo_u;
fe.v = clo_v;
}
if (x.u)
return this->valueFE(ItgPoint(x.u[0],x.u[1]));
return this->valueFE(fe);
// Use with caution, very slow!
Go::Point pt(x.x,x.y,x.z), clopt(3);
double clo_u, clo_v, dist;
#pragma omp critical
surf->closestPoint(pt, clo_u, clo_v, clopt, dist, 1.0e-5);
return this->valueFE(ItgPoint(clo_u,clo_v));
}
@ -160,7 +151,7 @@ bool SplineField2D::valueGrid (RealArray& val, const int* npe) const
}
bool SplineField2D::gradFE (const FiniteElement& fe, Vector& grad) const
bool SplineField2D::gradFE (const ItgPoint& x, Vector& grad) const
{
if (!basis) return false;
if (!surf) return false;
@ -168,7 +159,7 @@ bool SplineField2D::gradFE (const FiniteElement& fe, Vector& grad) const
// Evaluate the basis functions at the given point
Go::BasisDerivsSf spline;
#pragma omp critical
surf->computeBasis(fe.u,fe.v,spline);
surf->computeBasis(x.u,x.v,spline);
const int uorder = surf->order_u();
const int vorder = surf->order_v();
@ -189,14 +180,15 @@ bool SplineField2D::gradFE (const FiniteElement& fe, Vector& grad) const
Matrix Xnod(nsd,ip.size()), Jac;
for (size_t i = 0; i < ip.size(); i++)
Xnod.fillColumn(1+i,&(*surf->coefs_begin())+surf->dimension()*ip[i]);
utl::Jacobian(Jac,dNdX,Xnod,dNdu);
if (!utl::Jacobian(Jac,dNdX,Xnod,dNdu))
return false; // Singular Jacobian
// 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
#pragma omp critical
basis->computeBasis(fe.u,fe.v,spline);
basis->computeBasis(x.u,x.v,spline);
const size_t nbf = basis->order_u()*basis->order_v();
dNdu.resize(nbf,2);
@ -219,7 +211,7 @@ bool SplineField2D::gradFE (const FiniteElement& fe, Vector& grad) const
}
bool SplineField2D::hessianFE (const FiniteElement& fe, Matrix& H) const
bool SplineField2D::hessianFE (const ItgPoint& x, Matrix& H) const
{
if (!basis) return false;
if (!surf) return false;
@ -229,7 +221,7 @@ bool SplineField2D::hessianFE (const FiniteElement& fe, Matrix& H) const
IntVec ip;
if (surf == basis) {
#pragma omp critical
surf->computeBasis(fe.u,fe.v,spline2);
surf->computeBasis(x.u,x.v,spline2);
const size_t nen = surf->order_u()*surf->order_v();
d2Ndu2.resize(nen,2,2);
@ -246,7 +238,7 @@ bool SplineField2D::hessianFE (const FiniteElement& fe, Matrix& H) const
else {
// Mixed formulation, the solution uses a different basis than the geometry
#pragma omp critical
basis->computeBasis(fe.u,fe.v,spline2);
basis->computeBasis(x.u,x.v,spline2);
const size_t nbf = basis->order_u()*basis->order_v();
d2Ndu2.resize(nbf,2,2);

15
src/ASM/SplineField2D.h
View File

@ -37,7 +37,8 @@ public:
//! \param[in] patch The spline patch on which the field is to be defined
//! \param[in] v Array of control point field values
//! \param[in] basis Basis to use from patch
//! \param[in] cmp Component to use from source field. Pass 0 to use vector as-is.
//! \param[in] cmp Component to use from source field.
//! Pass 0 to use vector as-is.
//! \param[in] name Name of spline field
SplineField2D(const ASMs2D* patch, const RealArray& v,
char basis = 1, char cmp = 1, const char* name = nullptr);
@ -52,8 +53,8 @@ public:
virtual double valueNode(size_t node) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
virtual double valueFE(const FiniteElement& fe) const;
//! \param[in] x Local coordinate of evaluation point
virtual double valueFE(const ItgPoint& x) const;
//! \brief Computes the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point
@ -65,14 +66,14 @@ public:
virtual bool valueGrid(RealArray& val, const int* npe) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[in] x Local coordinate of evaluation point
//! \param[out] grad Gradient of solution in a given local coordinate
virtual bool gradFE(const FiniteElement& fe, Vector& grad) const;
virtual bool gradFE(const ItgPoint& x, Vector& grad) const;
//! \brief Computes the hessian for a given local coordinate.
//! \param[in] fe Finite element quantities
//! \param[in] x Local coordinate of evaluation point
//! \param[out] H Hessian of solution in a given local coordinate
virtual bool hessianFE(const FiniteElement& fe, Matrix& H) const;
virtual bool hessianFE(const ItgPoint& x, Matrix& H) const;
protected:
const Go::SplineSurface* basis; //!< Spline basis description

57
src/ASM/SplineField3D.C
View File

@ -15,7 +15,7 @@
#include "SplineField3D.h"
#include "ASMs3D.h"
#include "FiniteElement.h"
#include "ItgPoint.h"
#include "CoordinateMapping.h"
#include "Utilities.h"
#include "Vec3.h"
@ -61,14 +61,14 @@ double SplineField3D::valueNode (size_t node) const
}
double SplineField3D::valueFE (const FiniteElement& fe) const
double SplineField3D::valueFE (const ItgPoint& x) const
{
if (!basis) return false;
// Evaluate the basis functions at the given point
Go::BasisPts spline;
#pragma omp critical
basis->computeBasis(fe.u,fe.v,fe.w,spline);
basis->computeBasis(x.u,x.v,x.w,spline);
// Evaluate the solution field at the given point
IntVec ip;
@ -84,28 +84,16 @@ double SplineField3D::valueFE (const FiniteElement& fe) const
double SplineField3D::valueCoor (const Vec4& x) const
{
FiniteElement fe;
if (x.u) {
fe.u = x.u[0];
fe.v = x.u[1];
fe.w = x.u[2];
}
else {
// use with caution
Go::Point pt(3), clopt(3);
pt[0] = x[0];
pt[1] = x[1];
pt[2] = x[2];
double clo_u, clo_v, clo_w, dist;
if (x.u)
return this->valueFE(ItgPoint(x.u[0],x.u[1],x.u[2]));
// Use with caution, very slow!
Go::Point pt(x.x,x.y,x.z), clopt(3);
double clo_u, clo_v, clo_w, dist;
#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, 1.0e-5);
fe.u = clo_u;
fe.v = clo_v;
fe.w = clo_w;
}
return this->valueFE(fe);
return this->valueFE(ItgPoint(clo_u,clo_v,clo_w));
}
@ -143,13 +131,13 @@ bool SplineField3D::valueGrid (RealArray& val, const int* npe) const
// Evaluate the field in the visualization points
val.reserve(gpar[0].size()*gpar[1].size()*gpar[2].size());
for (size_t k = 0; k < gpar[2].size(); k++)
for (size_t j = 0; j < gpar[1].size(); j++)
for (size_t i = 0; i < gpar[0].size(); i++)
for (double w : gpar[2])
for (double v : gpar[1])
for (double u : gpar[0])
{
Go::BasisPts spline;
#pragma omp critical
basis->computeBasis(gpar[0][i],gpar[1][j],gpar[2][k],spline);
basis->computeBasis(u,v,w,spline);
IntVec ip;
ASMs3D::scatterInd(basis->numCoefs(0),basis->numCoefs(1),
@ -165,7 +153,7 @@ bool SplineField3D::valueGrid (RealArray& val, const int* npe) const
}
bool SplineField3D::gradFE (const FiniteElement& fe, Vector& grad) const
bool SplineField3D::gradFE (const ItgPoint& x, Vector& grad) const
{
if (!basis) return false;
if (!vol) return false;
@ -173,7 +161,7 @@ bool SplineField3D::gradFE (const FiniteElement& fe, Vector& grad) const
// Evaluate the basis functions at the given point
Go::BasisDerivs spline;
#pragma omp critical
vol->computeBasis(fe.u,fe.v,fe.w,spline);
vol->computeBasis(x.u,x.v,x.w,spline);
const int uorder = vol->order(0);
const int vorder = vol->order(1);
@ -196,14 +184,15 @@ bool SplineField3D::gradFE (const FiniteElement& fe, Vector& grad) const
Matrix Xnod(nsd,ip.size()), Jac;
for (size_t i = 0; i < ip.size(); i++)
Xnod.fillColumn(1+i,&(*vol->coefs_begin())+vol->dimension()*ip[i]);
utl::Jacobian(Jac,dNdX,Xnod,dNdu);
if (!utl::Jacobian(Jac,dNdX,Xnod,dNdu))
return false; // Singular Jacobian
// 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
#pragma omp critical
basis->computeBasis(fe.u,fe.v,fe.w,spline);
basis->computeBasis(x.u,x.v,x.w,spline);
const size_t nbf = basis->order(0)*basis->order(1)*basis->order(2);
dNdu.resize(nbf,3);
@ -227,7 +216,7 @@ bool SplineField3D::gradFE (const FiniteElement& fe, Vector& grad) const
}
bool SplineField3D::hessianFE (const FiniteElement& fe, Matrix& H) const
bool SplineField3D::hessianFE (const ItgPoint& x, Matrix& H) const
{
if (!basis) return false;
if (!vol) return false;
@ -237,7 +226,7 @@ bool SplineField3D::hessianFE (const FiniteElement& fe, Matrix& H) const
IntVec ip;
if (vol == basis) {
#pragma omp critical
vol->computeBasis(fe.u,fe.v,fe.w,spline2);
vol->computeBasis(x.u,x.v,x.w,spline2);
const size_t nen = vol->order(0)*vol->order(1)*vol->order(2);
d2Ndu2.resize(nen,3,3);
@ -257,7 +246,7 @@ bool SplineField3D::hessianFE (const FiniteElement& fe, Matrix& H) const
else {
// Mixed formulation, the solution uses a different basis than the geometry
#pragma omp critical
basis->computeBasis(fe.u,fe.v,fe.w,spline2);
basis->computeBasis(x.u,x.v,x.w,spline2);
const size_t nbf = basis->order(0)*basis->order(1)*basis->order(2);
d2Ndu2.resize(nbf,3,3);

17
src/ASM/SplineField3D.h
View File

@ -37,7 +37,8 @@ public:
//! \param[in] patch The spline patch on which the field is to be defined
//! \param[in] v Array of control point field values
//! \param[in] basis Basis to use from patch
//! \param[in] cmp Component to use from source field. Pass 0 to use vector as-is.
//! \param[in] cmp Component to use from source field.
//! Pass 0 to use the vector as-is.
//! \param[in] name Name of spline field
SplineField3D(const ASMs3D* patch, const RealArray& v,
char basis = 1, char cmp = 1, const char* name = nullptr);
@ -52,11 +53,11 @@ public:
virtual double valueNode(size_t node) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
virtual double valueFE(const FiniteElement& fe) const;
//! \param[in] x Local coordinate of evaluation point
virtual double valueFE(const ItgPoint& x) const;
//! \brief Computes the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point
//! \param[in] x Global/physical coordinate of evaluation point
virtual double valueCoor(const Vec4& x) const;
//! \brief Computes the value at a grid of visualization points.
@ -65,14 +66,14 @@ public:
virtual bool valueGrid(RealArray& val, const int* npe) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[in] x Local coordinate of evaluation point
//! \param[out] grad Gradient of solution in a given local coordinate
virtual bool gradFE(const FiniteElement& fe, Vector& grad) const;
virtual bool gradFE(const ItgPoint& x, Vector& grad) const;
//! \brief Computes the hessian for a given local coordinate.
//! \param[in] fe Finite element quantities
//! \param[in] x Local coordinate of evaluation point
//! \param[out] H Hessian of solution in a given local coordinate
virtual bool hessianFE(const FiniteElement& fe, Matrix& H) const;
virtual bool hessianFE(const ItgPoint& x, Matrix& H) const;
protected:
const Go::SplineVolume* basis; //!< Spline basis description

25
src/ASM/SplineFields1D.C
View File

@ -14,7 +14,7 @@
#include "GoTools/geometry/SplineCurve.h"
#include "SplineFields1D.h"
#include "FiniteElement.h"
#include "ItgPoint.h"
#include "CoordinateMapping.h"
#include "Utilities.h"
#include <numeric>
@ -30,13 +30,13 @@ SplineFields1D::SplineFields1D (const Go::SplineCurve* crv,
}
bool SplineFields1D::valueFE (const FiniteElement& fe, Vector& vals) const
bool SplineFields1D::valueFE (const ItgPoint& x, Vector& vals) const
{
if (!curv) return false;
// Find nodal indices of element containing given point
int first, last;
curv->basis().coefsAffectingParam(fe.u,first,last);
curv->basis().coefsAffectingParam(x.u,first,last);
if (first == last)
{
// We are at an interpolatory point (with C^0 continuity)
@ -54,7 +54,7 @@ bool SplineFields1D::valueFE (const FiniteElement& fe, Vector& vals) const
// Evaluate the basis functions at the given point
RealArray basisVal, basisDerivs;
#pragma omp critical
curv->computeBasis(fe.u,basisVal,basisDerivs);
curv->computeBasis(x.u,basisVal,basisDerivs);
// Evaluate the field at the given point.
// Notice we don't just do a matrix-vector multiplication here,
@ -70,20 +70,20 @@ bool SplineFields1D::valueFE (const FiniteElement& fe, Vector& vals) const
}
bool SplineFields1D::gradFE (const FiniteElement& fe, Matrix& grad) const
bool SplineFields1D::gradFE (const ItgPoint& x, Matrix& grad) const
{
if (!curv) return false;
// Evaluate the basis functions at the given point
RealArray basisVal, basisDerivs;
#pragma omp critical
curv->computeBasis(fe.u,basisVal,basisDerivs);
curv->computeBasis(x.u,basisVal,basisDerivs);
Matrix Jac, dNdX, dNdu(basisDerivs.size(),1);
dNdu.fillColumn(1,basisDerivs);
// Find nodal indices of element containing the given point
int first, last;
curv->basis().coefsAffectingParam(fe.u,first,last);
curv->basis().coefsAffectingParam(x.u,first,last);
std::vector<int> ip(last-first+1);
std::iota(ip.begin(),ip.end(),first);
@ -100,12 +100,13 @@ bool SplineFields1D::gradFE (const FiniteElement& fe, Matrix& grad) const
// Evaluate the field gradient at the given point
utl::gather(ip,nf,values,Velm);
if (!utl::Jacobian(Jac,dNdX,Xelm,dNdu))
return false;
return false; // Singular Jacobian
return !grad.multiply(Velm,dNdX).empty(); // grad = Xelm * dNdX
}
bool SplineFields1D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
bool SplineFields1D::hessianFE (const ItgPoint& x, Matrix3D& H) const
{
if (!curv) return false;
@ -116,7 +117,7 @@ bool SplineFields1D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
// Can this be fixed in GoTools?
Go::SplineCurve* ccrv = const_cast<Go::SplineCurve*>(curv);
#pragma omp critical
ccrv->computeBasis(fe.u,basisVal,basisDerivs1,basisDerivs2);
ccrv->computeBasis(x.u,basisVal,basisDerivs1,basisDerivs2);
Matrix Jac, dNdX, dNdu(basisDerivs1.size(),1);
dNdu.fillColumn(1,basisDerivs1);
Matrix3D Hess, d2NdX2, d2Ndu2(basisDerivs2.size(),1,1);
@ -124,8 +125,8 @@ bool SplineFields1D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
// Find nodal indices of element containing the given point
int first, last;
curv->basis().coefsAffectingParam(fe.u,first,last);
std::vector<int> ip(last-first+1);;
curv->basis().coefsAffectingParam(x.u,first,last);
std::vector<int> ip(last-first+1);
std::iota(ip.begin(),ip.end(),first);
// Find nodal coordinates of element containing the given point

12
src/ASM/SplineFields1D.h
View File

@ -46,19 +46,19 @@ public:
//==============================
//! \brief Computes the value for a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[in] x Local coordinate of evaluation point
//! \param[out] vals Values at local point in given element
virtual bool valueFE(const FiniteElement& fe, Vector& vals) const;
virtual bool valueFE(const ItgPoint& x, Vector& vals) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[in] x Local coordinate of evaluation point
//! \param[out] grad Gradient at local point in given element
virtual bool gradFE(const FiniteElement& fe, Matrix& grad) const;
virtual bool gradFE(const ItgPoint& x, Matrix& grad) const;
//! \brief Computes the hessian for a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[in] x Local coordinate of evaluation point
//! \param[out] H Hessian at local point in given element
virtual bool hessianFE(const FiniteElement& fe, Matrix3D& H) const;
virtual bool hessianFE(const ItgPoint& x, Matrix3D& H) const;
protected:
const Go::SplineCurve* curv; //!< Spline geometry description

51
src/ASM/SplineFields2D.C
View File

@ -15,7 +15,7 @@
#include "SplineFields2D.h"
#include "ASMs2D.h"
#include "FiniteElement.h"
#include "ItgPoint.h"
#include "CoordinateMapping.h"
#include "Utilities.h"
#include "Vec3.h"
@ -79,14 +79,14 @@ bool SplineFields2D::valueNode (size_t node, Vector& vals) const
}
bool SplineFields2D::valueFE (const FiniteElement& fe, Vector& vals) const
bool SplineFields2D::valueFE (const ItgPoint& x, Vector& vals) const
{
if (!basis) return false;
// Evaluate the basis functions at the given point
Go::BasisPtsSf spline;
#pragma omp critical
basis->computeBasis(fe.u,fe.v,spline);
basis->computeBasis(x.u,x.v,spline);
// Evaluate the solution field at the given point
std::vector<int> ip;
@ -104,30 +104,20 @@ bool SplineFields2D::valueFE (const FiniteElement& fe, Vector& vals) const
bool SplineFields2D::valueCoor (const Vec4& x, Vector& vals) const
{
FiniteElement fe;
if (x.u) {
fe.u = x.u[0];
fe.v = x.u[1];
}
else {
// use with caution
Go::Point pt(3), clopt(3);
pt[0] = x[0];
pt[1] = x[1];
pt[2] = x[2];
double clo_u, clo_v, dist;
if (x.u)
return this->valueFE(ItgPoint(x.u[0],x.u[1]),vals);
// Use with caution, very slow!
Go::Point pt(x.x,x.y,x.z), clopt(3);
double clo_u, clo_v, dist;
#pragma omp critical
surf->closestPoint(pt, clo_u, clo_v, clopt, dist, 1e-5);
surf->closestPoint(pt, clo_u, clo_v, clopt, dist, 1.0e-5);
fe.u = clo_u;
fe.v = clo_v;
}
return this->valueFE(fe, vals);
return this->valueFE(ItgPoint(clo_u,clo_v),vals);
}
bool SplineFields2D::gradFE (const FiniteElement& fe, Matrix& grad) const
bool SplineFields2D::gradFE (const ItgPoint& x, Matrix& grad) const
{
if (!basis) return false;
if (!surf) return false;
@ -135,7 +125,7 @@ bool SplineFields2D::gradFE (const FiniteElement& fe, Matrix& grad) const
// Evaluate the basis functions at the given point
Go::BasisDerivsSf spline;
#pragma omp critical
surf->computeBasis(fe.u,fe.v,spline);
surf->computeBasis(x.u,x.v,spline);
const int uorder = surf->order_u();
const int vorder = surf->order_v();
@ -156,14 +146,15 @@ bool SplineFields2D::gradFE (const FiniteElement& fe, Matrix& grad) const
Matrix Xnod(nsd,ip.size()), Jac;
for (size_t i = 0; i < ip.size(); i++)
Xnod.fillColumn(1+i,&(*surf->coefs_begin())+surf->dimension()*ip[i]);
utl::Jacobian(Jac,dNdX,Xnod,dNdu);
if (!utl::Jacobian(Jac,dNdX,Xnod,dNdu))
return false; // Singular Jacobian
// 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
#pragma omp critical
basis->computeBasis(fe.u,fe.v,spline);
basis->computeBasis(x.u,x.v,spline);
const size_t nbf = basis->order_u()*basis->order_v();
dNdu.resize(nbf,2);
@ -181,13 +172,11 @@ bool SplineFields2D::gradFE (const FiniteElement& fe, Matrix& grad) const
}
utl::gather(ip,nf,values,Xnod);
grad.multiply(Xnod,dNdX); // grad = Xnod * dNdX
return true;
return !grad.multiply(Xnod,dNdX).empty(); // grad = Xnod * dNdX
}
bool SplineFields2D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
bool SplineFields2D::hessianFE (const ItgPoint& x, Matrix3D& H) const
{
if (!basis) return false;
if (!surf) return false;
@ -198,7 +187,7 @@ bool SplineFields2D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
if (surf == basis) {
#pragma omp critical
surf->computeBasis(fe.u,fe.v,spline2);
surf->computeBasis(x.u,x.v,spline2);
const size_t nen = surf->order_u()*surf->order_v();
d2Ndu2.resize(nen,2,2);
@ -215,7 +204,7 @@ bool SplineFields2D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
else {
// Mixed formulation, the solution uses a different basis than the geometry
#pragma omp critical
basis->computeBasis(fe.u,fe.v,spline2);
basis->computeBasis(x.u,x.v,spline2);
const size_t nbf = basis->order_u()*basis->order_v();
d2Ndu2.resize(nbf,2,2);

12
src/ASM/SplineFields2D.h
View File

@ -62,9 +62,9 @@ public:
virtual bool valueNode(size_t node, Vector& vals) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[in] x Local coordinate of evaluation point
//! \param[out] vals Values in local point in given element
virtual bool valueFE(const FiniteElement& fe, Vector& vals) const;
virtual bool valueFE(const ItgPoint& x, Vector& vals) const;
//! \brief Computes the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point
@ -72,14 +72,14 @@ public:
virtual bool valueCoor(const Vec4& x, Vector& vals) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[in] x Local coordinate of evaluation point
//! \param[out] grad Gradient of solution in a given local coordinate
virtual bool gradFE(const FiniteElement& fe, Matrix& grad) const;
virtual bool gradFE(const ItgPoint& x, Matrix& grad) const;
//! \brief Computes the hessian for a given local coordinate.
//! \param[in] fe Finite element quantities
//! \param[in] x Local coordinate of evaluation point
//! \param[out] H Hessian of solution in a given local coordinate
virtual bool hessianFE(const FiniteElement& fe, Matrix3D& H) const;
virtual bool hessianFE(const ItgPoint& x, Matrix3D& H) const;
protected:
const Go::SplineSurface* basis; //!< Spline basis description

19
src/ASM/SplineFields2Dmx.C
View File

@ -11,14 +11,14 @@
//!
//==============================================================================
#include "GoTools/geometry/SplineSurface.h"
#include "SplineFields2Dmx.h"
#include "ASMs2Dmx.h"
#include "FiniteElement.h"
#include "ItgPoint.h"
#include "CoordinateMapping.h"
#include "Utilities.h"
#include "GoTools/geometry/SplineSurface.h"
SplineFields2Dmx::SplineFields2Dmx (const ASMs2Dmx* patch,
const RealArray& v, char basis,
@ -48,7 +48,7 @@ bool SplineFields2Dmx::valueNode (size_t node, Vector& vals) const
}
bool SplineFields2Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
bool SplineFields2Dmx::valueFE (const ItgPoint& x, Vector& vals) const
{
if (!surf) return false;
@ -61,7 +61,7 @@ bool SplineFields2Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
Go::SplineSurface* basis = surf->getBasis(b);
Go::BasisPtsSf spline;
#pragma omp critical
basis->computeBasis(fe.u,fe.v,spline);
basis->computeBasis(x.u,x.v,spline);
// Evaluate the solution field at the given point
std::vector<int> ip;
@ -81,7 +81,7 @@ bool SplineFields2Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
}
bool SplineFields2Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
bool SplineFields2Dmx::gradFE (const ItgPoint& x, Matrix& grad) const
{
if (!surf) return false;
@ -89,7 +89,7 @@ bool SplineFields2Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
Go::BasisDerivsSf spline;
const Go::SplineSurface* gsurf = surf->getBasis(ASMmxBase::geoBasis);
#pragma omp critical
gsurf->computeBasis(fe.u,fe.v,spline);
gsurf->computeBasis(x.u,x.v,spline);
const int uorder = gsurf->order_u();
const int vorder = gsurf->order_v();
@ -110,7 +110,8 @@ bool SplineFields2Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
Matrix Xnod(surf->getNoSpaceDim(),ip.size()), Jac;
for (size_t i = 0; i < ip.size(); i++)
Xnod.fillColumn(1+i,&(*gsurf->coefs_begin())+gsurf->dimension()*ip[i]);
utl::Jacobian(Jac,dNdX,Xnod,dNdu);
if (!utl::Jacobian(Jac,dNdX,Xnod,dNdu))
return false; // Singular Jacobian
// Evaluate the gradient of the solution field at the given point
auto vit = values.begin();
@ -118,7 +119,7 @@ bool SplineFields2Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
for (int b : bases) {
const Go::SplineSurface* basis = surf->getBasis(b);
#pragma omp critical
basis->computeBasis(fe.u,fe.v,spline);
basis->computeBasis(x.u,x.v,spline);
const size_t nbf = basis->order_u()*basis->order_v();
dNdu.resize(nbf,2);

16
src/ASM/SplineFields2Dmx.h
View File

@ -27,8 +27,8 @@ namespace Go {
/*!
\brief Class for mixed spline-based finite element vector fields in 2D.
\details This class implements the methods required to evaluate a 2D
mixed spline vector field at a given point in parametrical or physical coordinates.
\details This class implements the methods required to evaluate a 2D mixed
spline vector field at a given point in parametrical or physical coordinates.
*/
class SplineFields2Dmx : public Fields
@ -40,7 +40,7 @@ public:
//! \param[in] basis Bases to use from patch
//! \param[in] name Name of spline field
SplineFields2Dmx(const ASMs2Dmx* patch, const RealArray& v,
char basis = 12, const char* name = nullptr);
char basis = 12, const char* name = nullptr);
//! \brief Empty destructor.
virtual ~SplineFields2Dmx() {}
@ -53,18 +53,18 @@ public:
bool valueNode(size_t node, Vector& vals) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[in] x Local coordinate of evaluation point
//! \param[out] vals Values in local point in given element
bool valueFE(const FiniteElement& fe, Vector& vals) const;
bool valueFE(const ItgPoint& x, Vector& vals) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[in] x Local coordinate of evaluation point
//! \param[out] grad Gradient of solution in a given local coordinate
bool gradFE(const FiniteElement& fe, Matrix& grad) const;
bool gradFE(const ItgPoint& x, Matrix& grad) const;
protected:
const ASMs2Dmx* surf; //!< Patch description
std::set<int> bases; //!< Bases to use
std::set<int> bases; //!< Bases to use
};
#endif

53
src/ASM/SplineFields3D.C
View File

@ -15,7 +15,7 @@
#include "SplineFields3D.h"
#include "ASMs3D.h"
#include "FiniteElement.h"
#include "ItgPoint.h"
#include "CoordinateMapping.h"
#include "Utilities.h"
#include "Vec3.h"
@ -80,14 +80,14 @@ bool SplineFields3D::valueNode (size_t node, Vector& vals) const
}
bool SplineFields3D::valueFE (const FiniteElement& fe, Vector& vals) const
bool SplineFields3D::valueFE (const ItgPoint& x, Vector& vals) const
{
if (!basis) return false;
// Evaluate the basis functions at the given point
Go::BasisPts spline;
#pragma omp critical
basis->computeBasis(fe.u,fe.v,fe.w,spline);
basis->computeBasis(x.u,x.v,x.w,spline);
// Evaluate the solution field at the given point
std::vector<int> ip;
@ -105,32 +105,20 @@ bool SplineFields3D::valueFE (const FiniteElement& fe, Vector& vals) const
bool SplineFields3D::valueCoor (const Vec4& x, Vector& vals) const
{
FiniteElement fe;
if (x.u) {
fe.u = x.u[0];
fe.v = x.u[1];
fe.w = x.u[2];
}
else {
// use with caution
Go::Point pt(3), clopt(3);
pt[0] = x[0];
pt[1] = x[1];
pt[2] = x[2];
double clo_u, clo_v, clo_w, dist;
if (x.u)
return this->valueFE(ItgPoint(x.u[0],x.u[1],x.u[2]),vals);
// Use with caution, very slow!
Go::Point pt(x.x,x.y,x.z), clopt(3);
double clo_u, clo_v, clo_w, dist;
#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, 1.0e-5);
fe.u = clo_u;
fe.v = clo_v;
fe.w = clo_w;
}
return this->valueFE(fe, vals);
return this->valueFE(ItgPoint(clo_u,clo_v,clo_w),vals);
}
bool SplineFields3D::gradFE (const FiniteElement& fe, Matrix& grad) const
bool SplineFields3D::gradFE (const ItgPoint& x, Matrix& grad) const
{
if (!basis) return false;
if (!vol) return false;
@ -138,7 +126,7 @@ bool SplineFields3D::gradFE (const FiniteElement& fe, Matrix& grad) const
// Evaluate the basis functions at the given point
Go::BasisDerivs spline;
#pragma omp critical
vol->computeBasis(fe.u,fe.v,fe.w,spline);
vol->computeBasis(x.u,x.v,x.w,spline);
const int uorder = vol->order(0);
const int vorder = vol->order(1);
@ -161,14 +149,15 @@ bool SplineFields3D::gradFE (const FiniteElement& fe, Matrix& grad) const
Matrix Xnod(nsd,ip.size()), Jac;
for (size_t i = 0; i < ip.size(); i++)
Xnod.fillColumn(1+i,&(*vol->coefs_begin())+vol->dimension()*ip[i]);
utl::Jacobian(Jac,dNdX,Xnod,dNdu);
if (!utl::Jacobian(Jac,dNdX,Xnod,dNdu))
return false; // Singular Jacobian
// 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
#pragma omp critical
basis->computeBasis(fe.u,fe.v,fe.w,spline);
basis->computeBasis(x.u,x.v,x.w,spline);
const size_t nbf = basis->order(0)*basis->order(1)*basis->order(2);
dNdu.resize(nbf,3);
@ -187,13 +176,11 @@ bool SplineFields3D::gradFE (const FiniteElement& fe, Matrix& grad) const
}
utl::gather(ip,nf,values,Xnod);
grad.multiply(Xnod,dNdX); // grad = Xnod * dNdX
return true;
return !grad.multiply(Xnod,dNdX).empty(); // grad = Xnod * dNdX
}
bool SplineFields3D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
bool SplineFields3D::hessianFE (const ItgPoint& x, Matrix3D& H) const
{
if (!basis) return false;
if (!vol) return false;
@ -204,7 +191,7 @@ bool SplineFields3D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
IntVec ip;
if (vol == basis) {
#pragma omp critical
vol->computeBasis(fe.u,fe.v,fe.w,spline2);
vol->computeBasis(x.u,x.v,x.w,spline2);
const size_t nen = vol->order(0)*vol->order(1)*vol->order(2);
d2Ndu2.resize(nen,3,3);
@ -224,7 +211,7 @@ bool SplineFields3D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
else {
// Mixed formulation, the solution uses a different basis than the geometry
#pragma omp critical
basis->computeBasis(fe.u,fe.v,fe.w,spline2);
basis->computeBasis(x.u,x.v,x.w,spline2);
const size_t nbf = basis->order(0)*basis->order(1)*basis->order(2);
d2Ndu2.resize(nbf,3,3);

12
src/ASM/SplineFields3D.h
View File

@ -62,9 +62,9 @@ public:
virtual bool valueNode(size_t node, Vector& vals) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[in] x Local coordinate of evaluation point
//! \param[out] vals Values in local point in given element
virtual bool valueFE(const FiniteElement& fe, Vector& vals) const;
virtual bool valueFE(const ItgPoint& x, Vector& vals) const;
//! \brief Computes the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point
@ -72,14 +72,14 @@ public:
virtual bool valueCoor(const Vec4& x, Vector& vals) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[in] x Local coordinate of evaluation point
//! \param[out] grad Gradient of solution in a given local coordinate
virtual bool gradFE(const FiniteElement& fe, Matrix& grad) const;
virtual bool gradFE(const ItgPoint& x, Matrix& grad) const;
//! \brief Computes the hessian for a given local coordinate.
//! \param[in] fe Finite element quantities
//! \param[in] x Local coordinate of evaluation point
//! \param[out] H Hessian of solution in a given local coordinate
virtual bool hessianFE(const FiniteElement& fe, Matrix3D& H) const;
virtual bool hessianFE(const ItgPoint& x, Matrix3D& H) const;
protected:
const Go::SplineVolume* basis; //!< Spline basis description

23
src/ASM/SplineFields3Dmx.C
View File

@ -11,14 +11,14 @@
//!
//==============================================================================
#include "GoTools/trivariate/SplineVolume.h"
#include "SplineFields3Dmx.h"
#include "ASMs3Dmx.h"
#include "FiniteElement.h"
#include "ItgPoint.h"
#include "CoordinateMapping.h"
#include "Utilities.h"
#include "GoTools/trivariate/SplineVolume.h"
SplineFields3Dmx::SplineFields3Dmx (const ASMs3Dmx* patch,
const RealArray& v, char basis,
@ -48,7 +48,7 @@ bool SplineFields3Dmx::valueNode (size_t node, Vector& vals) const
}
bool SplineFields3Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
bool SplineFields3Dmx::valueFE (const ItgPoint& x, Vector& vals) const
{
if (!svol) return false;
@ -61,10 +61,10 @@ bool SplineFields3Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
Go::SplineVolume* basis = svol->getBasis(b);
Go::BasisPts spline;
#pragma omp critical
basis->computeBasis(fe.u,fe.v,fe.w,spline);
basis->computeBasis(x.u,x.v,x.w,spline);
// Evaluate the solution field at the given point
std::vector<int> ip;
IntVec ip;
ASMs3D::scatterInd(basis->numCoefs(0),basis->numCoefs(1),
basis->numCoefs(2),
basis->order(0),basis->order(1),basis->order(2),
@ -82,7 +82,7 @@ bool SplineFields3Dmx::valueFE (const FiniteElement& fe, Vector& vals) const
}
bool SplineFields3Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
bool SplineFields3Dmx::gradFE (const ItgPoint& x, Matrix& grad) const
{
if (!svol) return false;
@ -90,7 +90,7 @@ bool SplineFields3Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
Go::BasisDerivs spline;
const Go::SplineVolume* gvol = svol->getBasis(ASMmxBase::geoBasis);
#pragma omp critical
gvol->computeBasis(fe.u,fe.v,fe.w,spline);
gvol->computeBasis(x.u,x.v,x.w,spline);
const int uorder = gvol->order(0);
const int vorder = gvol->order(1);
@ -105,7 +105,7 @@ bool SplineFields3Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
dNdu(n,3) = spline.basisDerivs_w[n-1];
}
std::vector<int> ip;
IntVec ip;
ASMs3D::scatterInd(gvol->numCoefs(0),gvol->numCoefs(1),gvol->numCoefs(2),
uorder,vorder,worder,spline.left_idx,ip);
@ -113,7 +113,8 @@ bool SplineFields3Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
Matrix Xnod(svol->getNoSpaceDim(),ip.size()), Jac;
for (size_t i = 0; i < ip.size(); i++)
Xnod.fillColumn(1+i,&(*gvol->coefs_begin())+gvol->dimension()*ip[i]);
utl::Jacobian(Jac,dNdX,Xnod,dNdu);
if (!utl::Jacobian(Jac,dNdX,Xnod,dNdu))
return false; // Singular Jacobian
// Evaluate the gradient of the solution field at the given point
auto vit = values.begin();
@ -121,7 +122,7 @@ bool SplineFields3Dmx::gradFE (const FiniteElement& fe, Matrix& grad) const
for (int b : bases) {
const Go::SplineVolume* basis = svol->getBasis(b);
#pragma omp critical
basis->computeBasis(fe.u,fe.v,fe.w,spline);
basis->computeBasis(x.u,x.v,x.w,spline);
const size_t nbf = basis->order(0)*basis->order(1)*basis->order(2);
dNdu.resize(nbf,3);

16
src/ASM/SplineFields3Dmx.h
View File

@ -27,8 +27,8 @@ namespace Go {
/*!
\brief Class for mixed spline-based finite element vector fields in 3D.
\details This class implements the methods required to evaluate a 3D
mixed spline vector field at a given point in parametrical or physical coordinates.
\details This class implements the methods required to evaluate a 3D mixed
spline vector field at a given point in parametrical or physical coordinates.
*/
class SplineFields3Dmx : public Fields
@ -40,7 +40,7 @@ public:
//! \param[in] basis Bases to use from patch
//! \param[in] name Name of spline field
SplineFields3Dmx(const ASMs3Dmx* patch, const RealArray& v,
char basis = 12, const char* name = nullptr);
char basis = 12, const char* name = nullptr);
//! \brief Empty destructor.
virtual ~SplineFields3Dmx() {}
@ -53,18 +53,18 @@ public:
bool valueNode(size_t node, Vector& vals) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[in] x Local coordinate of evaluation point
//! \param[out] vals Values in local point in given element
bool valueFE(const FiniteElement& fe, Vector& vals) const;
bool valueFE(const ItgPoint& x, Vector& vals) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[in] x Local coordinate of evaluation point
//! \param[out] grad Gradient of solution in a given local coordinate
bool gradFE(const FiniteElement& fe, Matrix& grad) const;
bool gradFE(const ItgPoint& x, Matrix& grad) const;
protected:
const ASMs3Dmx* svol; //!< Patch description
std::set<int> bases; //!< Bases to use
std::set<int> bases; //!< Bases to use
};
#endif