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

Big revision of the spline field classes to avoid dangerous copy-pasting of internal GoTools code. Added a static create method in Field/Fields scope such that the fluid solvers can refer to fields independently of the discretization method.

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1327 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo 2011-12-01 13:30:01 +00:00 committed by Knut Morten Okstad
parent 7d67391641
commit 7532977a13
26 changed files with 681 additions and 339 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/ASM/ASMs2D.h
View File

@ -76,8 +76,10 @@ public:
//! \brief Empty destructor.
virtual ~ASMs2D() {}
//! \brief Returns the spline surface representing this patch.
//! \brief Returns the spline surface representing the geometry of this patch.
Go::SplineSurface* getSurface() const { return surf; }
//! \brief Returns the spline surface representing the basis of this patch.
virtual Go::SplineSurface* getBasis() const { return surf; }
// Methods for model generation

2
src/ASM/ASMs2DLag.h
View File

@ -138,6 +138,8 @@ protected:
//! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
//! in one element
virtual bool getElementCoordinates(Matrix& X, int iel) const;
public:
//! \brief Returns a matrix with all nodal coordinates within the patch.
//! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
//! in the patch

3
src/ASM/ASMs2Dmx.h
View File

@ -39,6 +39,9 @@ public:
//! \brief Empty destructor.
virtual ~ASMs2Dmx() {}
//! \brief Returns the spline surface representing the basis of this patch.
virtual Go::SplineSurface* getBasis() const { return basis1; }
//! \brief Returns the spline surface representing this patch.
Go::SplineSurface* getSurface() const { return surf; }

2
src/ASM/ASMs3D.h
View File

@ -98,6 +98,8 @@ public:
//! \brief Returns the spline volume representing this patch.
Go::SplineVolume* getVolume() const { return svol; }
//! \brief Returns the spline volume representing the basis of this patch.
virtual Go::SplineVolume* getBasis() const { return svol; }
// Methods for model generation

2
src/ASM/ASMs3DLag.h
View File

@ -148,6 +148,8 @@ protected:
//! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
//! in one element
virtual bool getElementCoordinates(Matrix& X, int iel) const;
public:
//! \brief Returns a matrix with all nodal coordinates within the patch.
//! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
//! in the patch

3
src/ASM/ASMs3Dmx.h
View File

@ -42,6 +42,9 @@ public:
//! \brief Empty destructor.
virtual ~ASMs3Dmx() {}
//! \brief Returns the spline volume representing the basis of this patch.
virtual Go::SplineVolume* getBasis() const { return basis1; }
// Methods for model generation
// ============================

37
src/ASM/Field.C Normal file
View File

@ -0,0 +1,37 @@
// $Id$
//==============================================================================
//!
//! \file Field.C
//!
//! \date Mar 28 2011
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Base class for scalar fields.
//!
//==============================================================================
#include "SplineField2D.h"
#include "SplineField3D.h"
#include "LagrangeField2D.h"
#include "LagrangeField3D.h"
#include "ASMs2DLag.h"
#include "ASMs3DLag.h"
Field* Field::create (const ASMbase* pch, const RealArray& v, const char* name)
{
const ASMs2DLag* pl2 = dynamic_cast<const ASMs2DLag*>(pch);
if (pl2) return new LagrangeField2D(pl2,v,name);
const ASMs2D* ps2 = dynamic_cast<const ASMs2D*>(pch);
if (ps2) return new SplineField2D(ps2,v,name);
const ASMs3DLag* pl3 = dynamic_cast<const ASMs3DLag*>(pch);
if (pl3) return new LagrangeField3D(pl3,v,name);
const ASMs3D* ps3 = dynamic_cast<const ASMs3D*>(pch);
if (ps3) return new SplineField3D(ps3,v,name);
return NULL;
}

34
src/ASM/Field.h
View File

@ -17,6 +17,7 @@
#include "MatVec.h"
#include <string>
class ASMbase;
class FiniteElement;
class Vec3;
@ -36,7 +37,8 @@ protected:
//! \brief The constructor sets the number of space dimensions and field name.
//! \param[in] n Number of space dimensions (1, 2 or 3)
//! \param[in] name Name of field
Field(unsigned char n, char* name = 0) : nsd(n) { if (name) fname = name; }
Field(unsigned char n, const char* name = 0) : nsd(n), nelm(0), nno(0)
{ if (name) fname = name; }
public:
//! \brief Empty destructor.
@ -46,27 +48,27 @@ public:
unsigned char getNoSpaceDim() const { return nsd; }
//! \brief Returns number of elements.
int getNoElm() const { return nelm; }
size_t getNoElm() const { return nelm; }
//! \brief Returns number of control points.
int getNoNodes() const { return nno; }
//! \brief Returns number of nodal/control points.
size_t getNoNodes() const { return nno; }
//! \brief Returns name of field.
//! \brief Returns the name of field.
const char* getFieldName() const { return fname.c_str(); }
//! \brief Sets the name of the field.
void setFieldName(const char* name) { fname = name; }
//! \brief Creates a dynamically allocated field object.
//! \param[in] pch The spline patch on which the field is to be defined on
//! \param[in] v Array of nodal/control point field values
//! \param[in] name Name of field
static Field* create(const ASMbase* pch, const RealArray& v,
const char* name = NULL);
//! \brief Initializes the field values.
void fill(const Vector& vec) { values = vec; }
// Methods to compute field values
//================================
// Methods to evaluate the field
//==============================
//! \brief Computes the value in a given node/control point.
//! \param[in] node Node number
virtual double valueNode(int node) const = 0;
virtual double valueNode(size_t node) const = 0;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
@ -88,8 +90,8 @@ public:
protected:
unsigned char nsd; //!< Number of space dimensions
int nelm; //!< Number of elements/knot-spans
int nno; //!< Number of nodes/control points
size_t nelm; //!< Number of elements/knot-spans
size_t nno; //!< Number of nodes/control points
std::string fname; //!< Name of the field
Vector values; //!< Field values
};

38
src/ASM/Fields.C Normal file
View File

@ -0,0 +1,38 @@
// $Id$
//==============================================================================
//!
//! \file Fields.C
//!
//! \date Mar 28 2011
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Base class for vector fields.
//!
//==============================================================================
#include "SplineFields2D.h"
#include "SplineFields3D.h"
#include "LagrangeFields2D.h"
#include "LagrangeFields3D.h"
#include "ASMs2DLag.h"
#include "ASMs3DLag.h"
Fields* Fields::create (const ASMbase* pch, const RealArray& v,
const char* name)
{
const ASMs2DLag* pl2 = dynamic_cast<const ASMs2DLag*>(pch);
if (pl2) return new LagrangeFields2D(pl2,v,name);
const ASMs2D* ps2 = dynamic_cast<const ASMs2D*>(pch);
if (ps2) return new SplineFields2D(ps2,v,name);
const ASMs3DLag* pl3 = dynamic_cast<const ASMs3DLag*>(pch);
if (pl3) return new LagrangeFields3D(pl3,v,name);
const ASMs3D* ps3 = dynamic_cast<const ASMs3D*>(pch);
if (ps3) return new SplineFields3D(ps3,v,name);
return NULL;
}

19
src/ASM/Fields.h
View File

@ -17,6 +17,7 @@
#include "MatVec.h"
#include <string>
class ASMbase;
class FiniteElement;
class Vec3;
@ -58,20 +59,20 @@ public:
//! \brief Returns name of field.
const char* getFieldName() const { return fname.c_str(); }
//! \brief Sets the name of the field.
void setFieldName(const char* name) { fname = name; }
//! \brief Creates a dynamically allocated field object.
//! \param[in] pch The spline patch on which the field is to be defined on
//! \param[in] v Array of nodal/control point field values
//! \param[in] name Name of field
static Fields* create(const ASMbase* pch, const RealArray& v,
const char* name = NULL);
//! \brief Initializes the field values.
void fill(const Vector& vec) { values = vec; nf = values.size()/nno; }
// Methods to compute field values
//================================
// Methods to evaluate the field
//==============================
//! \brief Computes the value in a given node/control point.
//! \param[in] node Node number
//! \param[out] vals Node values
virtual bool valueNode(int node, Vector& vals) const = 0;
virtual bool valueNode(size_t node, Vector& vals) const = 0;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition

18
src/ASM/LagrangeField2D.C
View File

@ -12,24 +12,32 @@
//==============================================================================
#include "LagrangeField2D.h"
#include "ASMs2DLag.h"
#include "FiniteElement.h"
#include "Lagrange.h"
#include "CoordinateMapping.h"
#include "Vec3.h"
LagrangeField2D::LagrangeField2D (const Matrix& X, int nx, int ny,
int px, int py, char* name)
: Field(2,name), coord(X), n1(nx), n2(ny), p1(px), p2(py)
LagrangeField2D::LagrangeField2D (const ASMs2DLag* patch, const RealArray& v,
const char* name) : Field(2,name)
{
patch->getNodalCoordinates(coord);
patch->getSize(n1,n2);
patch->getOrder(p1,p2);
nno = n1*n2;
nelm = (n1-1)*(n2-1)/(p1*p2);
// Ensure the values array has compatible length, pad with zeros if necessary
values.resize(nno);
RealArray::const_iterator end = v.size() > nno ? v.begin()+nno : v.end();
std::copy(v.begin(),end,values.begin());
}
double LagrangeField2D::valueNode (int node) const
double LagrangeField2D::valueNode (size_t node) const
{
return values(node);
return node > 0 && node <= nno ? values(node) : 0.0;
}

19
src/ASM/LagrangeField2D.h
View File

@ -16,6 +16,8 @@
#include "Field.h"
class ASMs2DLag;
/*!
\brief Class for Lagrange-based finite element scalar fields in 2D.
@ -28,23 +30,20 @@ class LagrangeField2D : public Field
{
public:
//! \brief The constructor sets the number of space dimensions and fields.
//! \param[in] X Matrix of nodal coordinates
//! \param[in] nx Number of nodes in first parameter direction
//! \param[in] ny Number of nodes in second parameter direction
//! \param[in] px Element order in first parameter direction
//! \param[in] py Element order in second parameter direction
//! \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] name Name of field
LagrangeField2D(const Matrix& X,
int nx, int ny, int px, int py, char* name = NULL);
LagrangeField2D(const ASMs2DLag* patch, const RealArray& v,
const char* name = NULL);
//! \brief Empty destructor.
virtual ~LagrangeField2D() {}
// Methods to compute field values
//================================
// Methods to evaluate the field
//==============================
//! \brief Computes the value in a given node/control point.
//! \param[in] node Node number
double valueNode(int node) const;
double valueNode(size_t node) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition

18
src/ASM/LagrangeField3D.C
View File

@ -12,24 +12,32 @@
//==============================================================================
#include "LagrangeField3D.h"
#include "ASMs3DLag.h"
#include "FiniteElement.h"
#include "Lagrange.h"
#include "CoordinateMapping.h"
#include "Vec3.h"
LagrangeField3D::LagrangeField3D (const Matrix& X, int nx, int ny, int nz,
int px, int py, int pz, char* name)
: Field(3,name), coord(X), n1(nx), n2(ny), n3(nz), p1(px), p2(py), p3(pz)
LagrangeField3D::LagrangeField3D (const ASMs3DLag* patch, const RealArray& v,
const char* name) : Field(3,name)
{
patch->getNodalCoordinates(coord);
patch->getSize(n1,n2,n3);
patch->getOrder(p1,p2,p3);
nno = n1*n2*n3;
nelm = (n1-1)*(n2-1)*(n3-1)/(p1*p2*p3);
// Ensure the values array has compatible length, pad with zeros if necessary
values.resize(nno);
RealArray::const_iterator end = v.size() > nno ? v.begin()+nno : v.end();
std::copy(v.begin(),end,values.begin());
}
double LagrangeField3D::valueNode (int node) const
double LagrangeField3D::valueNode (size_t node) const
{
return values(node);
return node > 0 && node <= nno ? values(node) : 0.0;
}

23
src/ASM/LagrangeField3D.h
View File

@ -16,6 +16,8 @@
#include "Field.h"
class ASMs3DLag;
/*!
\brief Class for Lagrange-based finite element scalar fields in 3D.
@ -28,31 +30,26 @@ class LagrangeField3D : public Field
{
public:
//! \brief The constructor sets the number of space dimensions and fields.
//! \param[in] X Matrix of nodal coordinates
//! \param[in] nx Number of nodes in first parameter direction
//! \param[in] ny Number of nodes in second parameter direction
//! \param[in] nz Number of nodes in third parameter direction
//! \param[in] px Element order in first parameter direction
//! \param[in] py Element order in second parameter direction
//! \param[in] pz Element order in third parameter direction
//! \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] name Name of field
LagrangeField3D(const Matrix& X, int nx, int ny, int nz,
int px, int py, int pz, char* name = NULL);
LagrangeField3D(const ASMs3DLag* patch, const RealArray& v,
const char* name = NULL);
//! \brief Empty destructor.
virtual ~LagrangeField3D() {}
// Methods to compute field values
//================================
// Methods to evaluate the field
//==============================
//! \brief Computes the value in a given node/control point.
//! \param[in] node Node number
double valueNode(int node) const;
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;
//! \brief Computed the value at a given global coordinate.
//! \brief Computes the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point
double valueCoor(const Vec3& x) const;

19
src/ASM/LagrangeFields2D.C
View File

@ -12,24 +12,33 @@
//==============================================================================
#include "LagrangeFields2D.h"
#include "ASMs2DLag.h"
#include "FiniteElement.h"
#include "Lagrange.h"
#include "CoordinateMapping.h"
#include "Vec3.h"
LagrangeFields2D::LagrangeFields2D (const Matrix& X, int nx, int ny,
int px, int py, char* name)
: Fields(2,name), coord(X), n1(nx), n2(ny), p1(px), p2(py)
LagrangeFields2D::LagrangeFields2D (const ASMs2DLag* patch, const RealArray& v,
const char* name) : Fields(2,name)
{
patch->getNodalCoordinates(coord);
patch->getSize(n1,n2);
patch->getOrder(p1,p2);
nno = n1*n2;
nelm = (n1-1)*(n2-1)/(p1*p2);
nf = v.size()/nno;
// Ensure the values array has compatible length, pad with zeros if necessary
values.resize(nf*nno);
RealArray::const_iterator end = v.size() > nf*nno ? v.begin()+nf*nno:v.end();
std::copy(v.begin(),end,values.begin());
}
bool LagrangeFields2D::valueNode (int node, Vector& vals) const
bool LagrangeFields2D::valueNode (size_t node, Vector& vals) const
{
if (node < 1 || (size_t)node > nno) return false;
if (node < 1 || node > nno) return false;
vals.resize(nf);
vals.fill(values.ptr()+nf*(node-1));

21
src/ASM/LagrangeFields2D.h
View File

@ -16,6 +16,8 @@
#include "Fields.h"
class ASMs2DLag;
/*!
\brief Class for Lagrange-based finite element vector fields in 2D.
@ -27,25 +29,22 @@
class LagrangeFields2D : public Fields
{
public:
//! \brief The constructor sets the field name.
//! \param[in] X Matrix of nodal coordinates
//! \param[in] nx Number of nodes in first parameter direction
//! \param[in] ny Number of nodes in second parameter direction
//! \param[in] px Element order in first parameter direction
//! \param[in] py Element order in second parameter direction
//! \brief The constructor sets the number of space dimensions and fields.
//! \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] name Name of field
LagrangeFields2D(const Matrix& X,
int nx, int ny, int px, int py, char* name = NULL);
LagrangeFields2D(const ASMs2DLag* patch, const RealArray& v,
const char* name = NULL);
//! \brief Empty destructor.
virtual ~LagrangeFields2D() {}
// Methods to compute field values
//================================
// Methods to evaluate the field
//==============================
//! \brief Computes the value in a given node/control point.
//! \param[in] node Node number
//! \param[out] vals Node values
bool valueNode(int node, Vector& vals) const;
bool valueNode(size_t node, Vector& vals) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition

19
src/ASM/LagrangeFields3D.C
View File

@ -12,24 +12,33 @@
//==============================================================================
#include "LagrangeFields3D.h"
#include "ASMs3DLag.h"
#include "FiniteElement.h"
#include "Lagrange.h"
#include "CoordinateMapping.h"
#include "Vec3.h"
LagrangeFields3D::LagrangeFields3D (const Matrix& X, int nx, int ny, int nz,
int px, int py, int pz, char* name)
: Fields(3,name), coord(X), n1(nx), n2(ny), n3(nz), p1(px), p2(py), p3(pz)
LagrangeFields3D::LagrangeFields3D (const ASMs3DLag* patch, const RealArray& v,
const char* name) : Fields(3,name)
{
patch->getNodalCoordinates(coord);
patch->getSize(n1,n2,n3);
patch->getOrder(p1,p2,p3);
nno = n1*n2*n3;
nelm = (n1-1)*(n2-1)*(n3-1)/(p1*p2*p3);
nf = v.size()/nno;
// Ensure the values array has compatible length, pad with zeros if necessary
values.resize(nf*nno);
RealArray::const_iterator end = v.size() > nf*nno ? v.begin()+nf*nno:v.end();
std::copy(v.begin(),end,values.begin());
}
bool LagrangeFields3D::valueNode (int node, Vector& vals) const
bool LagrangeFields3D::valueNode (size_t node, Vector& vals) const
{
if (node < 1 || (size_t)node > nno) return false;
if (node < 1 || node > nno) return false;
vals.resize(nf);
vals.fill(values.ptr()+nf*(node-1));

19
src/ASM/LagrangeFields3D.h
View File

@ -16,6 +16,8 @@
#include "Fields.h"
class ASMs3DLag;
/*!
\brief Class for Lagrange-based finite element vector fields in 3D.
@ -27,17 +29,12 @@
class LagrangeFields3D : public Fields
{
public:
//! \brief The constructor sets the field name.
//! \param[in] X Matrix of nodal coordinates
//! \param[in] nx Number of nodes in first parameter direction
//! \param[in] ny Number of nodes in second parameter direction
//! \param[in] nz Number of nodes in third parameter direction
//! \param[in] px Element order in first parameter direction
//! \param[in] py Element order in second parameter direction
//! \param[in] pz Element order in third parameter direction
//! \brief The constructor sets the number of space dimensions and fields.
//! \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] name Name of field
LagrangeFields3D(const Matrix& X, int nx, int ny, int nz,
int px, int py, int pz, char* name = NULL);
LagrangeFields3D(const ASMs3DLag* patch, const RealArray& v,
const char* name = NULL);
//! \brief Empty destructor.
virtual ~LagrangeFields3D() {}
@ -47,7 +44,7 @@ public:
//! \brief Computes the value in a given node/control point.
//! \param[in] node Node number
//! \param[out] vals Node values
bool valueNode(int node, Vector& vals) const;
bool valueNode(size_t node, Vector& vals) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition

134
src/ASM/SplineField2D.C
View File

@ -12,15 +12,18 @@
//==============================================================================
#include "SplineField2D.h"
#include "ASMs2D.h"
#include "FiniteElement.h"
#include "CoordinateMapping.h"
#include "Utilities.h"
#include "Vec3.h"
#include "GoTools/geometry/SplineSurface.h"
#include "GoTools/geometry/SplineUtils.h"
SplineField2D::SplineField2D(Go::SplineSurface *bf, char* name)
: Field(2,name), basis(bf), surf(bf)
SplineField2D::SplineField2D (const ASMs2D* patch,
const RealArray& v, const char* name)
: Field(2,name), basis(patch->getBasis()), surf(patch->getSurface())
{
const int n1 = basis->numCoefs_u();
const int n2 = basis->numCoefs_v();
@ -29,44 +32,41 @@ SplineField2D::SplineField2D(Go::SplineSurface *bf, char* name)
const int p1 = basis->order_u();
const int p2 = basis->order_v();
nelm = (n1-p1+1)*(n2-p2+1);
// Ensure the values array has compatible length, pad with zeros if necessary
values.resize(nno);
RealArray::const_iterator end = v.size() > nno ? v.begin()+nno : v.end();
std::copy(v.begin(),end,values.begin());
}
SplineField2D::SplineField2D(Go::SplineSurface *bf,
Go::SplineSurface *geometry,
char* name)
: Field(2,name), basis(bf), surf(geometry)
double SplineField2D::valueNode (size_t node) const
{
const int n1 = basis->numCoefs_u();
const int n2 = basis->numCoefs_v();
nno = n1*n2;
const int p1 = basis->order_u();
const int p2 = basis->order_v();
nelm = (n1-p1+1)*(n2-p2+1);
return node > 0 && node <= nno ? values(node) : 0.0;
}
SplineField2D::~SplineField2D()
double SplineField2D::valueFE (const FiniteElement& fe) const
{
// Set basis and geometry pointers to NULL, should not be deallocated here
basis = surf = NULL;
if (!basis) return false;
// Evaluate the basis functions at the given point
Go::BasisPtsSf spline;
basis->computeBasis(fe.u,fe.v,spline);
// Evaluate the solution field at the given point
std::vector<int> ip;
ASMs2D::scatterInd(basis->numCoefs_u(),basis->numCoefs_v(),
basis->order_u(),basis->order_v(),
spline.left_idx,ip);
Vector Vnod;
utl::gather(ip,1,values,Vnod);
return Vnod.dot(spline.basisValues);
}
double SplineField2D::valueNode(int node) const
{
// Not implemented yet
return 0.0;
}
double SplineField2D::valueFE(const FiniteElement& fe) const
{
// Go::Point pt;
// basis->pointValue(pt,values,fe.u,fe.v);
// return pt[0];
/* Old procedure, way too complex...
The above code is more in line with how we do it in ASMs2D.
const int uorder = basis->order_u();
const int vorder = basis->order_v();
const int unum = basis->numCoefs_u();
@ -140,19 +140,82 @@ double SplineField2D::valueFE(const FiniteElement& fe) const
return val;
}
*/
double SplineField2D::valueCoor(const Vec3& x) const
double SplineField2D::valueCoor (const Vec3& x) const
{
// Not implemented yet
return 0.0;
}
bool SplineField2D::gradFE(const FiniteElement& fe, Vector& grad) const
bool SplineField2D::gradFE (const FiniteElement& fe, Vector& grad) const
{
if (!basis) return false;
if (!surf) return false;
// Evaluate the basis functions at the given point
/*
Go::BasisDerivsSf spline;
surf->computeBasis(fe.u,fe.v,spline);
TODO: The above is not available yet, the below is temporary workaround.
*/
std::vector<Go::BasisDerivsSf> tmpSpline(1);
surf->computeBasisGrid(RealArray(1,fe.u),RealArray(1,fe.v),tmpSpline);
Go::BasisDerivsSf& spline = tmpSpline.front();
const int uorder = surf->order_u();
const int vorder = surf->order_v();
const size_t nen = uorder*vorder;
Matrix dNdu(nen,2), dNdX;
for (size_t n = 1; n <= nen; n++)
{
dNdu(n,1) = spline.basisDerivs_u[n-1];
dNdu(n,2) = spline.basisDerivs_v[n-1];
}
std::vector<int> ip;
ASMs2D::scatterInd(surf->numCoefs_u(),surf->numCoefs_v(),
uorder,vorder,spline.left_idx,ip);
// Evaluate the Jacobian inverse
Matrix Xnod, Jac;
Vector Xctrl(&(*surf->coefs_begin()),surf->coefs_end()-surf->coefs_begin());
utl::gather(ip,surf->dimension(),Xctrl,Xnod);
utl::Jacobian(Jac,dNdX,Xnod,dNdu);
// Evaluate the gradient of the solution field at the given point
if (basis != surf)
{
// Mixed formulation, the solution uses a different basis than the geometry
/*
basis->computeBasis(fe.u,fe.v,spline);
TODO: The above is not available yet, the below is temporary workaround.
*/
basis->computeBasisGrid(RealArray(1,fe.u),RealArray(1,fe.v),tmpSpline);
Go::BasisDerivsSf& spline = tmpSpline.front();
const size_t nbf = basis->order_u()*basis->order_v();
dNdu.resize(nbf,2);
for (size_t n = 1; n <= nbf; n++)
{
dNdu(n,1) = spline.basisDerivs_u[n-1];
dNdu(n,2) = spline.basisDerivs_v[n-1];
}
dNdX.multiply(dNdu,Jac); // dNdX = dNdu * Jac
ASMs2D::scatterInd(basis->numCoefs_u(),basis->numCoefs_v(),
basis->order_u(),basis->order_v(),
spline.left_idx,ip);
}
Vector Vnod;
utl::gather(ip,1,values,Vnod);
return dNdX.multiply(Vnod,grad,true);
}
/*
// // Derivatives of solution in parametric space
// std::vector<Go::Point> pts(3);
// basis->pointValue(pts,values,fe.u,fe.v,1);
@ -303,9 +366,10 @@ bool SplineField2D::gradFE(const FiniteElement& fe, Vector& grad) const
return true;
}
*/
bool SplineField2D::gradCoor(const Vec3& x, Vector& grad) const
bool SplineField2D::gradCoor (const Vec3& x, Vector& grad) const
{
// Not implemented yet
return false;

32
src/ASM/SplineField2D.h
View File

@ -16,6 +16,8 @@
#include "Field.h"
class ASMs2D;
namespace Go {
class SplineSurface;
}
@ -24,40 +26,34 @@ namespace Go {
/*!
\brief Class for spline-based finite element scalar fields in 2D.
\details This class implements the functions required to evaluate a 2D
\details This class implements the methods required to evaluate a 2D
spline scalar field at a given point in parametrical or physical coordinates.
*/
class SplineField2D : public Field
{
public:
//! \brief The constructor sets the number of space dimensions and fields.
//! \param[in] bf Spline basis description
//! \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] name Name of spline field
SplineField2D(Go::SplineSurface *bf, char* name = NULL);
//! \brief The constructor sets the number of space dimensions and fields.
//! \param[in] bf Spline basis description
//! \param[in] geometry Spline geometry description
//! \param[in] name Name of spline field
SplineField2D(Go::SplineSurface *bf,
Go::SplineSurface *geometry,
char* name = NULL);
SplineField2D(const ASMs2D* patch, const RealArray& v,
const char* name = NULL);
//! \brief Empty destructor.
virtual ~SplineField2D();
virtual ~SplineField2D() {}
// Methods to compute field values
//================================
// Methods to evaluate the field
//==============================
//! \brief Computes the value in a given node/control point.
//! \param[in] node Node number
double valueNode(int node) const;
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;
//! \brief Computed the value at a given global coordinate.
//! \brief Computes the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point
double valueCoor(const Vec3& x) const;
@ -72,8 +68,8 @@ public:
bool gradCoor(const Vec3& x, Vector& grad) const;
protected:
Go::SplineSurface* basis; //!< Spline basis description
Go::SplineSurface* surf; //!< Spline geometry description
const Go::SplineSurface* basis; //!< Spline basis description
const Go::SplineSurface* surf; //!< Spline geometry description
};
#endif

144
src/ASM/SplineField3D.C
View File

@ -12,19 +12,18 @@
//==============================================================================
#include "SplineField3D.h"
#include "ASMs3D.h"
#include "FiniteElement.h"
#include "CoordinateMapping.h"
#include "Utilities.h"
#include "Vec3.h"
#include "GoTools/trivariate/SplineVolume.h"
#include "GoTools/geometry/SplineUtils.h"
namespace Go {
void volume_ratder(double const eder[],int idim,int ider,double gder[]);
}
SplineField3D::SplineField3D(Go::SplineVolume *bf, char* name)
: Field(3,name), basis(bf), vol(bf)
SplineField3D::SplineField3D (const ASMs3D* patch,
const RealArray& v, const char* name)
: Field(3,name), basis(patch->getBasis()), vol(patch->getVolume())
{
const int n1 = basis->numCoefs(0);
const int n2 = basis->numCoefs(1);
@ -35,46 +34,41 @@ SplineField3D::SplineField3D(Go::SplineVolume *bf, char* name)
const int p2 = basis->order(1);
const int p3 = basis->order(2);
nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1);
// Ensure the values array has compatible length, pad with zeros if necessary
values.resize(nno);
RealArray::const_iterator end = v.size() > nno ? v.begin()+nno : v.end();
std::copy(v.begin(),end,values.begin());
}
SplineField3D::SplineField3D(Go::SplineVolume *bf,
Go::SplineVolume *geometry,
char* name)
: Field(3,name), basis(bf), vol(geometry)
double SplineField3D::valueNode (size_t node) const
{
const int n1 = basis->numCoefs(0);
const int n2 = basis->numCoefs(1);
const int n3 = basis->numCoefs(2);
nno = n1*n2*n3;
const int p1 = basis->order(0);
const int p2 = basis->order(1);
const int p3 = basis->order(2);
nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1);
return node > 0 && node <= nno ? values(node) : 0.0;
}
SplineField3D::~SplineField3D()
double SplineField3D::valueFE (const FiniteElement& fe) const
{
// Set basis and geometry pointers to NULL, should not be deallocated here
basis = vol = NULL;
if (!basis) return false;
// Evaluate the basis functions at the given point
Go::BasisPts spline;
basis->computeBasis(fe.u,fe.v,fe.w,spline);
// Evaluate the solution field at the given point
std::vector<int> ip;
ASMs3D::scatterInd(basis->numCoefs(0),basis->numCoefs(1),basis->numCoefs(2),
basis->order(0),basis->order(1),basis->order(2),
spline.left_idx,ip);
Vector Vnod;
utl::gather(ip,1,values,Vnod);
return Vnod.dot(spline.basisValues);
}
double SplineField3D::valueNode(int node) const
{
// Not implemented yet
return 0.0;
}
double SplineField3D::valueFE(const FiniteElement& fe) const
{
// Go::Point pt;
// basis->pointValue(pt,values,fe.u,fe.v,fe.w);
// return pt[0];
/* Old procedure, way too complex...
The above code is more in line with how we do it in ASMs3D.
double val = 0.0;
const int uorder = basis->order(0);
@ -168,19 +162,85 @@ double SplineField3D::valueFE(const FiniteElement& fe) const
val = tempResult/w;
return val;
}
*/
double SplineField3D::valueCoor(const Vec3& x) const
double SplineField3D::valueCoor (const Vec3& x) const
{
// Not implemented yet
return 0.0;
}
bool SplineField3D::gradFE(const FiniteElement& fe, Vector& grad) const
bool SplineField3D::gradFE (const FiniteElement& fe, Vector& grad) const
{
if (!basis) return false;
if (!vol) return false;
// Evaluate the basis functions at the given point
/*
Go::BasisDerivs spline;
vol->computeBasis(fe.u,fe.v,fe.w,spline);
TODO: The above is not available yet, the below is temporary workaround.
*/
std::vector<Go::BasisDerivs> tmpSpline(1);
vol->computeBasisGrid(RealArray(1,fe.u),RealArray(1,fe.v),RealArray(1,fe.w),tmpSpline);
Go::BasisDerivs& spline = tmpSpline.front();
const int uorder = vol->order(0);
const int vorder = vol->order(1);
const int worder = vol->order(2);
const size_t nen = uorder*vorder*worder;
Matrix dNdu(nen,3), dNdX;
for (size_t n = 1; n <= nen; n++)
{
dNdu(n,1) = spline.basisDerivs_u[n-1];
dNdu(n,2) = spline.basisDerivs_v[n-1];
dNdu(n,3) = spline.basisDerivs_w[n-1];
}
std::vector<int> ip;
ASMs3D::scatterInd(vol->numCoefs(0),vol->numCoefs(1),vol->numCoefs(2),
uorder,vorder,worder,spline.left_idx,ip);
// Evaluate the Jacobian inverse
Matrix Xnod, Jac;
Vector Xctrl(&(*vol->coefs_begin()),vol->coefs_end()-vol->coefs_begin());
utl::gather(ip,vol->dimension(),Xctrl,Xnod);
utl::Jacobian(Jac,dNdX,Xnod,dNdu);
// Evaluate the gradient of the solution field at the given point
if (basis != vol)
{
// Mixed formulation, the solution uses a different basis than the geometry
/*
basis->computeBasis(fe.u,fe.v,fe.w,spline);
TODO: The above is not available yet, the below is temporary workaround.
*/
basis->computeBasisGrid(RealArray(1,fe.u),RealArray(1,fe.v),RealArray(1,fe.w),tmpSpline);
Go::BasisDerivs& spline = tmpSpline.front();
const size_t nbf = basis->order(0)*basis->order(1)*basis->order(2);
dNdu.resize(nbf,3);
for (size_t n = 1; n <= nbf; n++)
{
dNdu(n,1) = spline.basisDerivs_u[n-1];
dNdu(n,2) = spline.basisDerivs_v[n-1];
dNdu(n,3) = spline.basisDerivs_w[n-1];
}
dNdX.multiply(dNdu,Jac);
ASMs3D::scatterInd(basis->numCoefs(0),basis->numCoefs(1),basis->numCoefs(2),
basis->order(0),basis->order(1),basis->order(2),
spline.left_idx,ip);
}
Vector Vnod;
utl::gather(ip,1,values,Vnod);
return dNdX.multiply(Vnod,grad,true);
}
/*
// // Derivatives of solution in parametric space
// std::vector<Go::Point> pts;
// basis->pointValue(pts,values,fe.u,fe.v,fe.w,1);
@ -376,9 +436,9 @@ bool SplineField3D::gradFE(const FiniteElement& fe, Vector& grad) const
return true;
}
*/
bool SplineField3D::gradCoor(const Vec3& x, Vector& grad) const
bool SplineField3D::gradCoor (const Vec3& x, Vector& grad) const
{
// Not implemented yet
return false;

29
src/ASM/SplineField3D.h
View File

@ -16,6 +16,8 @@
#include "Field.h"
class ASMs3D;
namespace Go {
class SplineVolume;
}
@ -28,35 +30,30 @@ namespace Go {
spline scalar field at a given point in parametrical or physical coordinates.
*/
class SplineField3D : public Field
{
public:
//! \brief The constructor sets the number of space dimensions and fields.
//! \param[in] bf Spline basis description
//! \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] name Name of spline field
SplineField3D(Go::SplineVolume *bf, char* name = NULL);
//! \param[in] bf Spline basis description
//! \param[in] geometry Spline geometry description
//! \param[in] name Name of spline field
SplineField3D(Go::SplineVolume *bf,
Go::SplineVolume *geometry,
char* name = NULL);
SplineField3D(const ASMs3D* patch, const RealArray& v,
const char* name = NULL);
//! \brief Empty destructor.
virtual ~SplineField3D();
virtual ~SplineField3D() {}
// Methods to compute field values
//================================
// Methods to evaluate the field
//==============================
//! \brief Computes the value in a given node/control point.
//! \param[in] node Node number
double valueNode(int node) const;
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;
//! \brief Computed the value at a given global coordinate.
//! \brief Computes the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point
double valueCoor(const Vec3& x) const;
@ -71,8 +68,8 @@ public:
bool gradCoor(const Vec3& x, Vector& grad) const;
protected:
Go::SplineVolume* basis; //!< Spline basis description
Go::SplineVolume* vol; //!< Spline volume description
const Go::SplineVolume* basis; //!< Spline basis description
const Go::SplineVolume* vol; //!< Spline geometry description
};
#endif

150
src/ASM/SplineFields2D.C
View File

@ -12,18 +12,18 @@
//==============================================================================
#include "SplineFields2D.h"
#include "ASMs2D.h"
#include "FiniteElement.h"
#include "CoordinateMapping.h"
#include "Utilities.h"
#include "Vec3.h"
#include "GoTools/geometry/SplineSurface.h"
#include "GoTools/geometry/SplineUtils.h"
#ifndef __BORLANDC__
#include "boost/lambda/lambda.hpp"
#endif
SplineFields2D::SplineFields2D(Go::SplineSurface *bf, char* name)
: Fields(2,name), basis(bf), surf(bf)
SplineFields2D::SplineFields2D (const ASMs2D* patch,
const RealArray& v, const char* name)
: Fields(2,name), basis(patch->getBasis()), surf(patch->getSurface())
{
const int n1 = basis->numCoefs_u();
const int n2 = basis->numCoefs_v();
@ -33,53 +33,47 @@ SplineFields2D::SplineFields2D(Go::SplineSurface *bf, char* name)
const int p2 = basis->order_v();
nelm = (n1-p1+1)*(n2-p2+1);
// Number of fields set in fill
nf = 0;
// Ensure the values array has compatible length, pad with zeros if necessary
nf = v.size()/nno;
values.resize(nf*nno);
RealArray::const_iterator end = v.size() > nf*nno ? v.begin()+nf*nno:v.end();
std::copy(v.begin(),end,values.begin());
}
SplineFields2D::SplineFields2D(Go::SplineSurface *bf,
Go::SplineSurface *geometry,
char* name)
: Fields(2,name), basis(bf), surf(geometry)
bool SplineFields2D::valueNode (size_t node, Vector& vals) const
{
const int n1 = basis->numCoefs_u();
const int n2 = basis->numCoefs_v();
nno = n1*n2;
if (node < 1 || node > nno) return false;
const int p1 = basis->order_u();
const int p2 = basis->order_v();
nelm = (n1-p1+1)*(n2-p2+1);
// Number of fields set in fill
nf = 0;
vals.resize(nf);
vals.fill(values.ptr()+(node-1)*nf);
return true;
}
SplineFields2D::~SplineFields2D()
{
// Set geometry to NULL; should not be deallocated here
basis = surf = NULL;
}
bool SplineFields2D::valueNode(int node, Vector& vals) const
{
// Not implemented yet
return false;
}
bool SplineFields2D::valueFE(const FiniteElement& fe, Vector& vals) const
bool SplineFields2D::valueFE (const FiniteElement& fe, Vector& vals) const
{
if (!basis) return false;
// Go::Point pt;
// basis->pointValue(pt,values,fe.u,fe.v);
// vals.resize(pt.size());
// for (int i = 0;i < pt.size();i++)
// vals[i] = pt[i];
// Evaluate the basis functions at the given point
Go::BasisPtsSf spline;
basis->computeBasis(fe.u,fe.v,spline);
// Evaluate the solution field at the given point
std::vector<int> ip;
ASMs2D::scatterInd(basis->numCoefs_u(),basis->numCoefs_v(),
basis->order_u(),basis->order_v(),
spline.left_idx,ip);
Matrix Vnod;
utl::gather(ip,nf,values,Vnod);
Vnod.multiply(spline.basisValues,vals); // vals = Vnod * basisValues
return true;
}
/* Old procedure, way too complex...
The above code is more in line with how we do it in ASMs2D.
const int uorder = basis->order_u();
const int vorder = basis->order_v();
const int unum = basis->numCoefs_u();
@ -175,19 +169,82 @@ bool SplineFields2D::valueFE(const FiniteElement& fe, Vector& vals) const
return true;
}
*/
bool SplineFields2D::valueCoor(const Vec3& x, Vector& vals) const
bool SplineFields2D::valueCoor (const Vec3& x, Vector& vals) const
{
// Not implemented yet
return false;
}
bool SplineFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const
bool SplineFields2D::gradFE (const FiniteElement& fe, Matrix& grad) const
{
if (!basis) return false;
if (!surf) return false;
// Evaluate the basis functions at the given point
/*
Go::BasisDerivsSf spline;
surf->computeBasis(fe.u,fe.v,spline);
TODO: The above is not available yet, the below is temporary workaround.
*/
std::vector<Go::BasisDerivsSf> tmpSpline(1);
surf->computeBasisGrid(RealArray(1,fe.u),RealArray(1,fe.v),tmpSpline);
Go::BasisDerivsSf& spline = tmpSpline.front();
const int uorder = surf->order_u();
const int vorder = surf->order_v();
const size_t nen = uorder*vorder;
Matrix dNdu(nen,2), dNdX;
for (size_t n = 1; n <= nen; n++)
{
dNdu(n,1) = spline.basisDerivs_u[n-1];
dNdu(n,2) = spline.basisDerivs_v[n-1];
}
std::vector<int> ip;
ASMs2D::scatterInd(surf->numCoefs_u(),surf->numCoefs_v(),
uorder,vorder,spline.left_idx,ip);
// Evaluate the Jacobian inverse
Matrix Xnod, Jac;
Vector Xctrl(&(*surf->coefs_begin()),surf->coefs_end()-surf->coefs_begin());
utl::gather(ip,surf->dimension(),Xctrl,Xnod);
utl::Jacobian(Jac,dNdX,Xnod,dNdu);
// Evaluate the gradient of the solution field at the given point
if (basis != surf)
{
// Mixed formulation, the solution uses a different basis than the geometry
/*
basis->computeBasis(fe.u,fe.v,spline);
TODO: The above is not available yet, the below is temporary workaround.
*/
basis->computeBasisGrid(RealArray(1,fe.u),RealArray(1,fe.v),tmpSpline);
Go::BasisDerivsSf& spline = tmpSpline.front();
const size_t nbf = basis->order_u()*basis->order_v();
dNdu.resize(nbf,2);
for (size_t n = 1; n <= nbf; n++)
{
dNdu(n,1) = spline.basisDerivs_u[n-1];
dNdu(n,2) = spline.basisDerivs_v[n-1];
}
dNdX.multiply(dNdu,Jac); // dNdX = dNdu * Jac
ASMs2D::scatterInd(basis->numCoefs_u(),basis->numCoefs_v(),
basis->order_u(),basis->order_v(),
spline.left_idx,ip);
}
utl::gather(ip,nf,values,Xnod);
grad.multiply(Xnod,dNdX); // grad = Xnod * dNdX
return true;
}
// // Derivatives of solution in parametric space
// std::vector<Go::Point> pts(3);
@ -215,7 +272,7 @@ bool SplineFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const
// // Gradient of solution in given parametric point
// // df/dX = df/dXi * Jac^-1 = df/dXi * [dX/dXi]^-1
// grad.multiply(gradXi,Jac);
/*
// Gradient of field
Matrix gradXi(nf,2);
@ -351,9 +408,10 @@ bool SplineFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const
return true;
}
*/
bool SplineFields2D::gradCoor(const Vec3& x, Matrix& grad) const
bool SplineFields2D::gradCoor (const Vec3& x, Matrix& grad) const
{
// Not implemented yet
return false;

36
src/ASM/SplineFields2D.h
View File

@ -16,6 +16,8 @@
#include "Fields.h"
class ASMs2D;
namespace Go {
class SplineSurface;
}
@ -24,44 +26,38 @@ namespace Go {
/*!
\brief Class for spline-based finite element vector fields in 2D.
\details This class implements the functions required to evaluate a 2D
\details This class implements the methods required to evaluate a 2D
spline vector field at a given point in parametrical or physical coordinates.
*/
class SplineFields2D : public Fields
{
public:
//! \brief The constructor sets the field name.
//! \param[in] bf Spline basis description
//! \brief The constructor sets the number of space dimensions and fields.
//! \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] name Name of spline field
SplineFields2D(Go::SplineSurface *bf, char* name = NULL);
//! \brief The constructor sets the field name.
//! \param[in] bf Spline basis description
//! \param[in] geometry Spline geometry description
//! \param[in] name Name of spline field
SplineFields2D(Go::SplineSurface *bf,
Go::SplineSurface *geoemtry,
char* name = NULL);
SplineFields2D(const ASMs2D* patch, const RealArray& v,
const char* name = NULL);
//! \brief Empty destructor.
virtual ~SplineFields2D();
virtual ~SplineFields2D() {}
// Methods to compute field values
//================================
// Methods to evaluate the field
//==============================
//! \brief Computes the value in a given node/control point.
//! \param[in] node Node number
//! \param[out] vals Node values
bool valueNode(int node, Vector& vals) const;
bool valueNode(size_t node, Vector& vals) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[out] vals Values in local point in given element
bool valueFE(const FiniteElement& fe, Vector& vals) const;
//! \brief Computed the value at a given global coordinate.
//! \brief Computes the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point
//! \param[in] vals Values in given physical coordinate
//! \param[out] vals Values in given physical coordinate
bool valueCoor(const Vec3& x, Vector& vals) const;
//! \brief Computes the gradient for a given local coordinate.
@ -75,8 +71,8 @@ public:
bool gradCoor(const Vec3& x, Matrix& grad) const;
protected:
Go::SplineSurface* basis; //!< Spline basis description
Go::SplineSurface* surf; //!< Spline geoemtry description
const Go::SplineSurface* basis; //!< Spline basis description
const Go::SplineSurface* surf; //!< Spline geometry description
};
#endif

159
src/ASM/SplineFields3D.C
View File

@ -12,20 +12,18 @@
//==============================================================================
#include "SplineFields3D.h"
#include "ASMs3D.h"
#include "FiniteElement.h"
#include "CoordinateMapping.h"
#include "Utilities.h"
#include "Vec3.h"
#include "GoTools/trivariate/SplineVolume.h"
#include "GoTools/geometry/SplineUtils.h"
#include "boost/lambda/lambda.hpp"
namespace Go {
void volume_ratder(double const eder[],int idim,int ider,double gder[]);
}
SplineFields3D::SplineFields3D(Go::SplineVolume *bf, char* name)
: Fields(3,name), basis(bf), vol(bf)
SplineFields3D::SplineFields3D (const ASMs3D* patch,
const RealArray& v, const char* name)
: Fields(3,name), basis(patch->getBasis()), vol(patch->getVolume())
{
const int n1 = basis->numCoefs(0);
const int n2 = basis->numCoefs(1);
@ -37,56 +35,48 @@ SplineFields3D::SplineFields3D(Go::SplineVolume *bf, char* name)
const int p3 = basis->order(2);
nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1);
// Number of fields set in fill
nf = 0;
// Ensure the values array has compatible length, pad with zeros if necessary
nf = v.size()/nno;
values.resize(nf*nno);
RealArray::const_iterator end = v.size() > nf*nno ? v.begin()+nf*nno:v.end();
std::copy(v.begin(),end,values.begin());
}
SplineFields3D::SplineFields3D(Go::SplineVolume *bf,
Go::SplineVolume *geometry,
char* name)
: Fields(3,name), basis(bf), vol(geometry)
bool SplineFields3D::valueNode (size_t node, Vector& vals) const
{
const int n1 = basis->numCoefs(0);
const int n2 = basis->numCoefs(1);
const int n3 = basis->numCoefs(2);
nno = n1*n2*n3;
if (node < 1 || node > nno) return false;
const int p1 = basis->order(0);
const int p2 = basis->order(1);
const int p3 = basis->order(2);
nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1);
// Number of fields set in fill
nf = 0;
vals.resize(nf);
vals.fill(values.ptr()+(node-1)*nf);
return true;
}
SplineFields3D::~SplineFields3D()
{
// Set geometry to NULL; should not be deallocated here
basis = NULL;
}
bool SplineFields3D::valueNode(int node, Vector& vals) const
{
// Not implemented yet
return false;
}
bool SplineFields3D::valueFE(const FiniteElement& fe, Vector& vals) const
bool SplineFields3D::valueFE (const FiniteElement& fe, Vector& vals) const
{
if (!basis) return false;
Go::Point pt;
// Evaluate the basis functions at the given point
Go::BasisPts spline;
basis->computeBasis(fe.u,fe.v,fe.w,spline);
// basis->pointValue(pt,values,fe.u,fe.v,fe.w);
// vals.resize(pt.size());
// for (int i = 0;i < pt.size();i++)
// vals[i] = pt[i];
// Evaluate the solution field at the given point
std::vector<int> ip;
ASMs3D::scatterInd(basis->numCoefs(0),basis->numCoefs(1),basis->numCoefs(2),
basis->order(0),basis->order(1),basis->order(2),
spline.left_idx,ip);
Matrix Vnod;
utl::gather(ip,nf,values,Vnod);
Vnod.multiply(spline.basisValues,vals); // vals = Vnod * basisValues
return true;
}
/* Old procedure, way too complex...
The above code is more in line with how we do it in ASMs3D.
Go::Point pt;
const int uorder = basis->order(0);
const int vorder = basis->order(1);
@ -213,10 +203,10 @@ bool SplineFields3D::valueFE(const FiniteElement& fe, Vector& vals) const
return true;
}
*/
bool SplineFields3D::valueCoor(const Vec3& x, Vector& vals) const
bool SplineFields3D::valueCoor (const Vec3& x, Vector& vals) const
{
// Not implemented yet
return false;
@ -226,7 +216,73 @@ bool SplineFields3D::valueCoor(const Vec3& x, Vector& vals) const
bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const
{
if (!basis) return false;
if (!vol) return false;
// Evaluate the basis functions at the given point
/*
Go::BasisDerivs spline;
vol->computeBasis(fe.u,fe.v,fe.w,spline);
TODO: The above is not available yet, the below is temporary workaround.
*/
std::vector<Go::BasisDerivs> tmpSpline(1);
vol->computeBasisGrid(RealArray(1,fe.u),RealArray(1,fe.v),RealArray(1,fe.w),tmpSpline);
Go::BasisDerivs& spline = tmpSpline.front();
const int uorder = vol->order(0);
const int vorder = vol->order(1);
const int worder = vol->order(2);
const size_t nen = uorder*vorder*worder;
Matrix dNdu(nen,3), dNdX;
for (size_t n = 1; n <= nen; n++)
{
dNdu(n,1) = spline.basisDerivs_u[n-1];
dNdu(n,2) = spline.basisDerivs_v[n-1];
dNdu(n,3) = spline.basisDerivs_w[n-1];
}
std::vector<int> ip;
ASMs3D::scatterInd(vol->numCoefs(0),vol->numCoefs(1),vol->numCoefs(2),
uorder,vorder,worder,spline.left_idx,ip);
// Evaluate the Jacobian inverse
Matrix Xnod, Jac;
Vector Xctrl(&(*vol->coefs_begin()),vol->coefs_end()-vol->coefs_begin());
utl::gather(ip,vol->dimension(),Xctrl,Xnod);
utl::Jacobian(Jac,dNdX,Xnod,dNdu);
// Evaluate the gradient of the solution field at the given point
if (basis != vol)
{
// Mixed formulation, the solution uses a different basis than the geometry
/*
basis->computeBasis(fe.u,fe.v,fe.w,spline);
TODO: The above is not available yet, the below is temporary workaround.
*/
basis->computeBasisGrid(RealArray(1,fe.u),RealArray(1,fe.v),RealArray(1,fe.w),tmpSpline);
Go::BasisDerivs& spline = tmpSpline.front();
const size_t nbf = basis->order(0)*basis->order(1)*basis->order(2);
dNdu.resize(nbf,3);
for (size_t n = 1; n <= nbf; n++)
{
dNdu(n,1) = spline.basisDerivs_u[n-1];
dNdu(n,2) = spline.basisDerivs_v[n-1];
dNdu(n,3) = spline.basisDerivs_w[n-1];
}
dNdX.multiply(dNdu,Jac); // dNdX = dNdu * Jac
ASMs3D::scatterInd(basis->numCoefs(0),basis->numCoefs(1),basis->numCoefs(2),
basis->order(0),basis->order(1),basis->order(2),
spline.left_idx,ip);
}
utl::gather(ip,nf,values,Xnod);
grad.multiply(Xnod,dNdX); // grad = Xnod * dNdX
return true;
}
// // Derivatives of solution in parametric space
// std::vector<Go::Point> pts;
// basis->pointValue(pts,values,fe.u,fe.v,fe.w,1);
@ -252,6 +308,7 @@ bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const
// // df/dX = df/dXi * Jac^-1 = df/dXi * [dX/dXi]^-1
// grad.multiply(gradXi,Jac);
/*
// Derivatives of solution in parametric space
std::vector<Go::Point> pts(4);
@ -438,9 +495,9 @@ bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const
return true;
}
*/
bool SplineFields3D::gradCoor(const Vec3& x, Matrix& grad) const
bool SplineFields3D::gradCoor (const Vec3& x, Matrix& grad) const
{
// Not implemented yet
return false;

36
src/ASM/SplineFields3D.h
View File

@ -16,6 +16,8 @@
#include "Fields.h"
class ASMs3D;
namespace Go {
class SplineVolume;
}
@ -24,44 +26,38 @@ namespace Go {
/*!
\brief Class for spline-based finite element vector fields in 3D.
\details This class implements the functions required to evaluate a 3D
\details This class implements the methods required to evaluate a 3D
spline vector field at a given point in parametrical or physical coordinates.
*/
class SplineFields3D : public Fields
{
public:
//! \brief The constructor sets the field name.
//! \param[in] bf Spline basis description
//! \brief The constructor sets the number of space dimensions and fields.
//! \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] name Name of spline field
SplineFields3D(Go::SplineVolume *bf, char* name = NULL);
//! \brief The constructor sets the field name.
//! \param[in] bf Spline basis description
//! \param[in] geometry Spline geometry description
//! \param[in] name Name of spline field
SplineFields3D(Go::SplineVolume *bf,
Go::SplineVolume *geometry,
char* name = NULL);
SplineFields3D(const ASMs3D* patch, const RealArray& v,
const char* name = NULL);
//! \brief Empty destructor.
virtual ~SplineFields3D();
virtual ~SplineFields3D() {}
// Methods to compute field values
//================================
// Methods to evaluate the field
//==============================
//! \brief Computes the value in a given node/control point.
//! \param[in] node Node number
//! \param[out] vals Node values
bool valueNode(int node, Vector& vals) const;
bool valueNode(size_t node, Vector& vals) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[out] vals Values in local point in given element
bool valueFE(const FiniteElement& fe, Vector& vals) const;
//! \brief Computed the value at a given global coordinate.
//! \brief Computes the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point
//! \param[in] vals Values in given physical coordinate
//! \param[out] vals Values in given physical coordinate
bool valueCoor(const Vec3& x, Vector& vals) const;
//! \brief Computes the gradient for a given local coordinate.
@ -75,8 +71,8 @@ public:
bool gradCoor(const Vec3& x, Matrix& grad) const;
protected:
Go::SplineVolume* basis; //!< Spline basis description
Go::SplineVolume* vol; //!< Spline geometry description
const Go::SplineVolume* basis; //!< Spline basis description
const Go::SplineVolume* vol; //!< Spline geometry description
};
#endif