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

Added support for mixed formulation

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1303 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
rho 2011-11-29 10:42:40 +00:00 committed by Knut Morten Okstad
parent c5f34c7fb0
commit 24b88f6a25
8 changed files with 278 additions and 176 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

85
src/ASM/SplineField2D.C
View File

@ -19,23 +19,38 @@
#include "GoTools/geometry/SplineUtils.h" #include "GoTools/geometry/SplineUtils.h"
SplineField2D::SplineField2D(Go::SplineSurface *geometry, char* name) SplineField2D::SplineField2D(Go::SplineSurface *bf, char* name)
: Field(2,name), surf(geometry) : Field(2,name), basis(bf), surf(bf)
{ {
const int n1 = surf->numCoefs_u(); const int n1 = basis->numCoefs_u();
const int n2 = surf->numCoefs_v(); const int n2 = basis->numCoefs_v();
nno = n1*n2; nno = n1*n2;
const int p1 = surf->order_u(); const int p1 = basis->order_u();
const int p2 = surf->order_v(); const int p2 = basis->order_v();
nelm = (n1-p1+1)*(n2-p2+1);
}
SplineField2D::SplineField2D(Go::SplineSurface *bf,
Go::SplineSurface *geometry,
char* name)
: Field(2,name), basis(bf), surf(geometry)
{
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); nelm = (n1-p1+1)*(n2-p2+1);
} }
SplineField2D::~SplineField2D() SplineField2D::~SplineField2D()
{ {
// Set geometry pointer to NULL, should not be deallocated here // Set basis and geometry pointers to NULL, should not be deallocated here
surf = NULL; basis = surf = NULL;
} }
@ -49,16 +64,16 @@ double SplineField2D::valueNode(int node) const
double SplineField2D::valueFE(const FiniteElement& fe) const double SplineField2D::valueFE(const FiniteElement& fe) const
{ {
// Go::Point pt; // Go::Point pt;
// surf->pointValue(pt,values,fe.u,fe.v); // basis->pointValue(pt,values,fe.u,fe.v);
// return pt[0]; // return pt[0];
const int uorder = surf->order_u(); const int uorder = basis->order_u();
const int vorder = surf->order_v(); const int vorder = basis->order_v();
const int unum = surf->numCoefs_u(); const int unum = basis->numCoefs_u();
//const int vnum = surf->numCoefs_v(); //const int vnum = basis->numCoefs_v();
const int dim = surf->dimension(); const int dim = basis->dimension();
const bool rational = surf->rational(); const bool rational = basis->rational();
const int kdim = rational ? dim + 1 : dim; const int kdim = rational ? dim + 1 : dim;
//const int ncomp = values.size()/(unum*vnum); //const int ncomp = values.size()/(unum*vnum);
@ -71,11 +86,11 @@ double SplineField2D::valueFE(const FiniteElement& fe) const
Bu.resize(uorder); Bu.resize(uorder);
Bv.resize(vorder); Bv.resize(vorder);
surf->basis_u().computeBasisValues(fe.u,Bu.begin()); basis->basis_u().computeBasisValues(fe.u,Bu.begin());
surf->basis_v().computeBasisValues(fe.v,Bv.begin()); basis->basis_v().computeBasisValues(fe.v,Bv.begin());
const int uleft = surf->basis_u().lastKnotInterval(); const int uleft = basis->basis_u().lastKnotInterval();
const int vleft = surf->basis_v().lastKnotInterval(); const int vleft = basis->basis_v().lastKnotInterval();
// compute the tensor product value // compute the tensor product value
const int val_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1)); const int val_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1));
@ -86,7 +101,7 @@ double SplineField2D::valueFE(const FiniteElement& fe) const
double tempw; double tempw;
double w = 0.0; double w = 0.0;
const int w_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1)) * kdim + dim; const int w_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1)) * kdim + dim;
register const double* w_ptr = &(surf->rcoefs_begin()[w_start_ix]); register const double* w_ptr = &(basis->rcoefs_begin()[w_start_ix]);
for (register double* bval_v_ptr = Bv.begin(); bval_v_ptr != Bv.end(); ++bval_v_ptr) { for (register double* bval_v_ptr = Bv.begin(); bval_v_ptr != Bv.end(); ++bval_v_ptr) {
register const double bval_v = *bval_v_ptr; register const double bval_v = *bval_v_ptr;
@ -136,11 +151,11 @@ double SplineField2D::valueCoor(const Vec3& x) const
bool SplineField2D::gradFE(const FiniteElement& fe, Vector& grad) const bool SplineField2D::gradFE(const FiniteElement& fe, Vector& grad) const
{ {
if (!surf) return false; if (!basis) return false;
// // Derivatives of solution in parametric space // // Derivatives of solution in parametric space
// std::vector<Go::Point> pts(3); // std::vector<Go::Point> pts(3);
// surf->pointValue(pts,values,fe.u,fe.v,1); // basis->pointValue(pts,values,fe.u,fe.v,1);
// // Gradient of field wrt parametric coordinates // // Gradient of field wrt parametric coordinates
// Vector gradXi(2); // Vector gradXi(2);
@ -149,7 +164,7 @@ bool SplineField2D::gradFE(const FiniteElement& fe, Vector& grad) const
// // Derivatives of coordinate mapping // // Derivatives of coordinate mapping
// std::vector<Go::Point> xpts(3); // std::vector<Go::Point> xpts(3);
// surf->point(xpts,fe.u,fe.v,1); // basis->point(xpts,fe.u,fe.v,1);
// // Jacobian matrix // // Jacobian matrix
// Matrix Jac(2,2); // Matrix Jac(2,2);
@ -167,9 +182,9 @@ bool SplineField2D::gradFE(const FiniteElement& fe, Vector& grad) const
// Gradient of field wrt parametric coordinates // Gradient of field wrt parametric coordinates
Vector gradXi(2); Vector gradXi(2);
const int uorder = surf->order_u(); const int uorder = basis->order_u();
const int vorder = surf->order_v(); const int vorder = basis->order_v();
const int unum = surf->numCoefs_u(); const int unum = basis->numCoefs_u();
const int derivs = 1; const int derivs = 1;
const int totpts = 3; const int totpts = 3;
@ -178,11 +193,11 @@ bool SplineField2D::gradFE(const FiniteElement& fe, Vector& grad) const
const double resolution = 1.0e-12; const double resolution = 1.0e-12;
// Dimension // Dimension
const int dim = surf->dimension(); const int dim = basis->dimension();
const bool rational = surf->rational(); const bool rational = basis->rational();
// Take care of the rational case // Take care of the rational case
const std::vector<double>::iterator co = rational ? surf->rcoefs_begin() : surf->coefs_begin(); const std::vector<double>::iterator co = rational ? basis->rcoefs_begin() : basis->coefs_begin();
int kdim = dim + (rational ? 1 : 0); int kdim = dim + (rational ? 1 : 0);
int ndim = 1 + (rational ? 1 : 0); int ndim = 1 + (rational ? 1 : 0);
@ -196,16 +211,16 @@ bool SplineField2D::gradFE(const FiniteElement& fe, Vector& grad) const
// Compute the basis values and get some data about the spline spaces // Compute the basis values and get some data about the spline spaces
if (u_from_right) if (u_from_right)
surf->basis_u().computeBasisValues(fe.u, &b0[0], derivs, resolution); basis->basis_u().computeBasisValues(fe.u, &b0[0], derivs, resolution);
else else
surf->basis_u().computeBasisValuesLeft(fe.u, &b0[0], derivs, resolution); basis->basis_u().computeBasisValuesLeft(fe.u, &b0[0], derivs, resolution);
int uleft = surf->basis_u().lastKnotInterval(); int uleft = basis->basis_u().lastKnotInterval();
if (v_from_right) if (v_from_right)
surf->basis_v().computeBasisValues(fe.v, &b1[0], derivs, resolution); basis->basis_v().computeBasisValues(fe.v, &b1[0], derivs, resolution);
else else
surf->basis_v().computeBasisValuesLeft(fe.v, &b1[0], derivs, resolution); basis->basis_v().computeBasisValuesLeft(fe.v, &b1[0], derivs, resolution);
int vleft = surf->basis_v().lastKnotInterval(); int vleft = basis->basis_v().lastKnotInterval();
// Compute the tensor product value // Compute the tensor product value
int coefind = uleft-uorder+1 + unum*(vleft-vorder+1); int coefind = uleft-uorder+1 + unum*(vleft-vorder+1);

14
src/ASM/SplineField2D.h
View File

@ -33,9 +33,16 @@ class SplineField2D : public Field
{ {
public: public:
//! \brief The constructor sets the number of space dimensions and fields. //! \brief The constructor sets the number of space dimensions and fields.
//! \param[in] geometry Spline surface geometry //! \param[in] bf Spline basis description
//! \param[in] name Name of spline field //! \param[in] name Name of spline field
SplineField2D(Go::SplineSurface *geometry, char* name = NULL); 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);
//! \brief Empty destructor. //! \brief Empty destructor.
virtual ~SplineField2D(); virtual ~SplineField2D();
@ -65,7 +72,8 @@ public:
bool gradCoor(const Vec3& x, Vector& grad) const; bool gradCoor(const Vec3& x, Vector& grad) const;
protected: protected:
Go::SplineSurface* surf; //!< Spline surface geometry description Go::SplineSurface* basis; //!< Spline basis description
Go::SplineSurface* surf; //!< Spline geometry description
}; };
#endif #endif

109
src/ASM/SplineField3D.C
View File

@ -23,25 +23,42 @@ namespace Go {
} }
SplineField3D::SplineField3D(Go::SplineVolume *geometry, char* name) SplineField3D::SplineField3D(Go::SplineVolume *bf, char* name)
: Field(3,name), vol(geometry) : Field(3,name), basis(bf), vol(bf)
{ {
const int n1 = vol->numCoefs(0); const int n1 = basis->numCoefs(0);
const int n2 = vol->numCoefs(1); const int n2 = basis->numCoefs(1);
const int n3 = vol->numCoefs(2); const int n3 = basis->numCoefs(2);
nno = n1*n2*n3; nno = n1*n2*n3;
const int p1 = vol->order(0); const int p1 = basis->order(0);
const int p2 = vol->order(1); const int p2 = basis->order(1);
const int p3 = vol->order(2); const int p3 = basis->order(2);
nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1);
}
SplineField3D::SplineField3D(Go::SplineVolume *bf,
Go::SplineVolume *geometry,
char* name)
: Field(3,name), basis(bf), vol(geometry)
{
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); nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1);
} }
SplineField3D::~SplineField3D() SplineField3D::~SplineField3D()
{ {
// Set geometry pointer to NULL, should not be deallocated here // Set basis and geometry pointers to NULL, should not be deallocated here
vol = NULL; basis = vol = NULL;
} }
@ -55,19 +72,19 @@ double SplineField3D::valueNode(int node) const
double SplineField3D::valueFE(const FiniteElement& fe) const double SplineField3D::valueFE(const FiniteElement& fe) const
{ {
// Go::Point pt; // Go::Point pt;
// vol->pointValue(pt,values,fe.u,fe.v,fe.w); // basis->pointValue(pt,values,fe.u,fe.v,fe.w);
// return pt[0]; // return pt[0];
double val = 0.0; double val = 0.0;
const int uorder = vol->order(0); const int uorder = basis->order(0);
const int vorder = vol->order(1); const int vorder = basis->order(1);
const int worder = vol->order(2); const int worder = basis->order(2);
const int unum = vol->numCoefs(0); const int unum = basis->numCoefs(0);
const int vnum = vol->numCoefs(1); const int vnum = basis->numCoefs(1);
const int dim = vol->dimension(); const int dim = basis->dimension();
const bool rational = vol->rational(); const bool rational = basis->rational();
int kdim = rational ? dim + 1 : dim; int kdim = rational ? dim + 1 : dim;
static Go::ScratchVect<double, 10> Bu(uorder); static Go::ScratchVect<double, 10> Bu(uorder);
@ -80,12 +97,12 @@ double SplineField3D::valueFE(const FiniteElement& fe) const
Bw.resize(worder); Bw.resize(worder);
// compute tbe basis values and get some data about the spline spaces // compute tbe basis values and get some data about the spline spaces
vol->basis(0).computeBasisValues(fe.u, Bu.begin()); basis->basis(0).computeBasisValues(fe.u, Bu.begin());
vol->basis(1).computeBasisValues(fe.v, Bv.begin()); basis->basis(1).computeBasisValues(fe.v, Bv.begin());
vol->basis(2).computeBasisValues(fe.w, Bw.begin()); basis->basis(2).computeBasisValues(fe.w, Bw.begin());
const int uleft = vol->basis(0).lastKnotInterval(); const int uleft = basis->basis(0).lastKnotInterval();
const int vleft = vol->basis(1).lastKnotInterval(); const int vleft = basis->basis(1).lastKnotInterval();
const int wleft = vol->basis(2).lastKnotInterval(); const int wleft = basis->basis(2).lastKnotInterval();
// compute the tensor product value // compute the tensor product value
const int val_start_ix = uleft - uorder + 1 + unum * (vleft - vorder + 1 + vnum * (wleft - worder + 1)); const int val_start_ix = uleft - uorder + 1 + unum * (vleft - vorder + 1 + vnum * (wleft - worder + 1));
@ -97,7 +114,7 @@ double SplineField3D::valueFE(const FiniteElement& fe) const
w = 0.0; w = 0.0;
const int w_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1 + vnum * (wleft - worder + 1))) * kdim + dim; const int w_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1 + vnum * (wleft - worder + 1))) * kdim + dim;
const double *w_ptr = &(vol->rcoefs_begin()[w_start_ix]); const double *w_ptr = &(basis->rcoefs_begin()[w_start_ix]);
for (double* bval_w_ptr = Bw.begin(); bval_w_ptr != Bw.end(); ++bval_w_ptr) { for (double* bval_w_ptr = Bw.begin(); bval_w_ptr != Bw.end(); ++bval_w_ptr) {
const double bval_w = *bval_w_ptr; const double bval_w = *bval_w_ptr;
@ -162,11 +179,11 @@ double SplineField3D::valueCoor(const Vec3& x) const
bool SplineField3D::gradFE(const FiniteElement& fe, Vector& grad) const bool SplineField3D::gradFE(const FiniteElement& fe, Vector& grad) const
{ {
if (!vol) return false; if (!basis) return false;
// // Derivatives of solution in parametric space // // Derivatives of solution in parametric space
// std::vector<Go::Point> pts; // std::vector<Go::Point> pts;
// vol->pointValue(pts,values,fe.u,fe.v,fe.w,1); // basis->pointValue(pts,values,fe.u,fe.v,fe.w,1);
// // Gradient of field wrt parametric coordinates // // Gradient of field wrt parametric coordinates
// Vector gradXi(3); // Vector gradXi(3);
@ -176,7 +193,7 @@ bool SplineField3D::gradFE(const FiniteElement& fe, Vector& grad) const
// // Derivatives of coordinate mapping // // Derivatives of coordinate mapping
// std::vector<Go::Point> xpts(3); // std::vector<Go::Point> xpts(3);
// vol->point(xpts,fe.u,fe.v,fe.w,1); // basis->point(xpts,fe.u,fe.v,fe.w,1);
// // Jacobian matrix // // Jacobian matrix
// Matrix Jac(3,3); // Matrix Jac(3,3);
@ -197,13 +214,13 @@ bool SplineField3D::gradFE(const FiniteElement& fe, Vector& grad) const
int derivs = 1; int derivs = 1;
int totpts = (derivs + 1)*(derivs + 2)*(derivs + 3)/6; int totpts = (derivs + 1)*(derivs + 2)*(derivs + 3)/6;
const int unum = vol->numCoefs(0); const int unum = basis->numCoefs(0);
const int vnum = vol->numCoefs(1); const int vnum = basis->numCoefs(1);
const int wnum = vol->numCoefs(2); const int wnum = basis->numCoefs(2);
const int uorder = vol->order(0); const int uorder = basis->order(0);
const int vorder = vol->order(1); const int vorder = basis->order(1);
const int worder = vol->order(2); const int worder = basis->order(2);
const bool u_from_right = true; const bool u_from_right = true;
const bool v_from_right = true; const bool v_from_right = true;
@ -212,11 +229,11 @@ bool SplineField3D::gradFE(const FiniteElement& fe, Vector& grad) const
const double resolution = 1.0e-12; const double resolution = 1.0e-12;
const int ncomp = values.size()/(unum*vnum*wnum); const int ncomp = values.size()/(unum*vnum*wnum);
const int dim = vol->dimension(); const int dim = basis->dimension();
const bool rational = vol->rational(); const bool rational = basis->rational();
// Take care of the rational case // Take care of the rational case
const std::vector<double>::iterator co = rational ? vol->rcoefs_begin() : vol->coefs_begin(); const std::vector<double>::iterator co = rational ? basis->rcoefs_begin() : basis->coefs_begin();
int kdim = dim + (rational ? 1 : 0); int kdim = dim + (rational ? 1 : 0);
int ndim = ncomp + (rational ? 1 : 0); int ndim = ncomp + (rational ? 1 : 0);
@ -232,23 +249,23 @@ bool SplineField3D::gradFE(const FiniteElement& fe, Vector& grad) const
// Compute the basis values and get some data about the spline spaces // Compute the basis values and get some data about the spline spaces
if (u_from_right) if (u_from_right)
vol->basis(0).computeBasisValues(fe.u, &b0[0], derivs, resolution); basis->basis(0).computeBasisValues(fe.u, &b0[0], derivs, resolution);
else else
vol->basis(1).computeBasisValuesLeft(fe.u, &b0[0], derivs, resolution); basis->basis(1).computeBasisValuesLeft(fe.u, &b0[0], derivs, resolution);
int uleft = vol->basis(0).lastKnotInterval(); int uleft = basis->basis(0).lastKnotInterval();
if (v_from_right) if (v_from_right)
vol->basis(1).computeBasisValues(fe.v, &b1[0], derivs, resolution); basis->basis(1).computeBasisValues(fe.v, &b1[0], derivs, resolution);
else else
vol->basis(1).computeBasisValuesLeft(fe.v, &b1[0], derivs, resolution); basis->basis(1).computeBasisValuesLeft(fe.v, &b1[0], derivs, resolution);
int vleft = vol->basis(1).lastKnotInterval(); int vleft = basis->basis(1).lastKnotInterval();
if (w_from_right) if (w_from_right)
vol->basis(2).computeBasisValues(fe.w, &b2[0], derivs, resolution); basis->basis(2).computeBasisValues(fe.w, &b2[0], derivs, resolution);
else else
vol->basis(2).computeBasisValuesLeft(fe.w, &b2[0], derivs, resolution); basis->basis(2).computeBasisValuesLeft(fe.w, &b2[0], derivs, resolution);
int wleft = vol->basis(2).lastKnotInterval(); int wleft = basis->basis(2).lastKnotInterval();
// Compute the tensor product value // Compute the tensor product value
int coefind = uleft-uorder+1 + unum*(vleft-vorder+1 + vnum*(wleft-worder+1)); int coefind = uleft-uorder+1 + unum*(vleft-vorder+1 + vnum*(wleft-worder+1));

13
src/ASM/SplineField3D.h
View File

@ -33,9 +33,15 @@ class SplineField3D : public Field
{ {
public: public:
//! \brief The constructor sets the number of space dimensions and fields. //! \brief The constructor sets the number of space dimensions and fields.
//! \param[in] geometry Spline volume geometry //! \param[in] bf Spline basis description
//! \param[in] name Name of spline field //! \param[in] name Name of spline field
SplineField3D(Go::SplineVolume *geometry, char* name = NULL); 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);
//! \brief Empty destructor. //! \brief Empty destructor.
virtual ~SplineField3D(); virtual ~SplineField3D();
@ -65,7 +71,8 @@ public:
bool gradCoor(const Vec3& x, Vector& grad) const; bool gradCoor(const Vec3& x, Vector& grad) const;
protected: protected:
Go::SplineVolume* vol; //!< Spline volume geometry description Go::SplineVolume* basis; //!< Spline basis description
Go::SplineVolume* vol; //!< Spline volume description
}; };
#endif #endif

90
src/ASM/SplineFields2D.C
View File

@ -22,15 +22,33 @@
#endif #endif
SplineFields2D::SplineFields2D(Go::SplineSurface *geometry, char* name) SplineFields2D::SplineFields2D(Go::SplineSurface *bf, char* name)
: Fields(2,name), surf(geometry) : Fields(2,name), basis(bf), surf(bf)
{ {
const int n1 = surf->numCoefs_u(); const int n1 = basis->numCoefs_u();
const int n2 = surf->numCoefs_v(); const int n2 = basis->numCoefs_v();
nno = n1*n2; nno = n1*n2;
const int p1 = surf->order_u(); const int p1 = basis->order_u();
const int p2 = surf->order_v(); const int p2 = basis->order_v();
nelm = (n1-p1+1)*(n2-p2+1);
// Number of fields set in fill
nf = 0;
}
SplineFields2D::SplineFields2D(Go::SplineSurface *bf,
Go::SplineSurface *geometry,
char* name)
: Fields(2,name), basis(bf), surf(geometry)
{
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); nelm = (n1-p1+1)*(n2-p2+1);
// Number of fields set in fill // Number of fields set in fill
@ -41,7 +59,7 @@ SplineFields2D::SplineFields2D(Go::SplineSurface *geometry, char* name)
SplineFields2D::~SplineFields2D() SplineFields2D::~SplineFields2D()
{ {
// Set geometry to NULL; should not be deallocated here // Set geometry to NULL; should not be deallocated here
surf = NULL; basis = surf = NULL;
} }
@ -54,21 +72,21 @@ bool SplineFields2D::valueNode(int node, Vector& vals) const
bool SplineFields2D::valueFE(const FiniteElement& fe, Vector& vals) const bool SplineFields2D::valueFE(const FiniteElement& fe, Vector& vals) const
{ {
if (!surf) return false; if (!basis) return false;
// Go::Point pt; // Go::Point pt;
// surf->pointValue(pt,values,fe.u,fe.v); // basis->pointValue(pt,values,fe.u,fe.v);
// vals.resize(pt.size()); // vals.resize(pt.size());
// for (int i = 0;i < pt.size();i++) // for (int i = 0;i < pt.size();i++)
// vals[i] = pt[i]; // vals[i] = pt[i];
const int uorder = surf->order_u(); const int uorder = basis->order_u();
const int vorder = surf->order_v(); const int vorder = basis->order_v();
const int unum = surf->numCoefs_u(); const int unum = basis->numCoefs_u();
const int vnum = surf->numCoefs_v(); const int vnum = basis->numCoefs_v();
const int dim = surf->dimension(); const int dim = basis->dimension();
const bool rational = surf->rational(); const bool rational = basis->rational();
int kdim = rational ? dim + 1 : dim; int kdim = rational ? dim + 1 : dim;
const int ncomp = values.size()/(unum*vnum); const int ncomp = values.size()/(unum*vnum);
@ -83,11 +101,11 @@ bool SplineFields2D::valueFE(const FiniteElement& fe, Vector& vals) const
Bv.resize(vorder); Bv.resize(vorder);
tempVal.resize(ncomp); tempVal.resize(ncomp);
surf->basis_u().computeBasisValues(fe.u,Bu.begin()); basis->basis_u().computeBasisValues(fe.u,Bu.begin());
surf->basis_v().computeBasisValues(fe.v,Bv.begin()); basis->basis_v().computeBasisValues(fe.v,Bv.begin());
const int uleft = surf->basis_u().lastKnotInterval(); const int uleft = basis->basis_u().lastKnotInterval();
const int vleft = surf->basis_v().lastKnotInterval(); const int vleft = basis->basis_v().lastKnotInterval();
//compute the tensor product value //compute the tensor product value
const int val_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1)) * ncomp; const int val_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1)) * ncomp;
@ -101,7 +119,7 @@ bool SplineFields2D::valueFE(const FiniteElement& fe, Vector& vals) const
double w = 0.0; double w = 0.0;
const int w_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1)) * kdim + dim; const int w_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1)) * kdim + dim;
register const double* w_ptr = &(surf->rcoefs_begin()[w_start_ix]); register const double* w_ptr = &(basis->rcoefs_begin()[w_start_ix]);
for (register double* bval_v_ptr = Bv.begin(); bval_v_ptr != Bv.end(); ++bval_v_ptr) { for (register double* bval_v_ptr = Bv.begin(); bval_v_ptr != Bv.end(); ++bval_v_ptr) {
register const double bval_v = *bval_v_ptr; register const double bval_v = *bval_v_ptr;
@ -169,11 +187,11 @@ bool SplineFields2D::valueCoor(const Vec3& x, Vector& vals) const
bool SplineFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const bool SplineFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const
{ {
if (!surf) return false; if (!basis) return false;
// // Derivatives of solution in parametric space // // Derivatives of solution in parametric space
// std::vector<Go::Point> pts(3); // std::vector<Go::Point> pts(3);
// surf->pointValue(pts,values,fe.u,fe.v,1); // basis->pointValue(pts,values,fe.u,fe.v,1);
// // Gradient of field wrt parametric coordinates // // Gradient of field wrt parametric coordinates
// Matrix gradXi(nf,2); // Matrix gradXi(nf,2);
@ -183,7 +201,7 @@ bool SplineFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const
// // Derivatives of coordinate mapping // // Derivatives of coordinate mapping
// std::vector<Go::Point> xpts(3); // std::vector<Go::Point> xpts(3);
// surf->point(xpts,fe.u,fe.v,1); // basis->point(xpts,fe.u,fe.v,1);
// // Jacobian matrix // // Jacobian matrix
// Matrix Jac(2,2); // Matrix Jac(2,2);
@ -201,10 +219,10 @@ bool SplineFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const
// Gradient of field // Gradient of field
Matrix gradXi(nf,2); Matrix gradXi(nf,2);
const int uorder = surf->order_u(); const int uorder = basis->order_u();
const int vorder = surf->order_v(); const int vorder = basis->order_v();
const int unum = surf->numCoefs_u(); const int unum = basis->numCoefs_u();
const int vnum = surf->numCoefs_v(); const int vnum = basis->numCoefs_v();
const int derivs = 1; const int derivs = 1;
const bool u_from_right = true; const bool u_from_right = true;
@ -215,11 +233,11 @@ bool SplineFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const
const int ncomp = values.size()/(unum*vnum); const int ncomp = values.size()/(unum*vnum);
// Dimension // Dimension
const int dim = surf->dimension(); const int dim = basis->dimension();
const bool rational = surf->rational(); const bool rational = basis->rational();
// Take care of the rational case // Take care of the rational case
const std::vector<double>::iterator co = rational ? surf->rcoefs_begin() : surf->coefs_begin(); const std::vector<double>::iterator co = rational ? basis->rcoefs_begin() : basis->coefs_begin();
int kdim = dim + (rational ? 1 : 0); int kdim = dim + (rational ? 1 : 0);
int ndim = ncomp + (rational ? 1 : 0); int ndim = ncomp + (rational ? 1 : 0);
@ -233,18 +251,18 @@ bool SplineFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const
// Compute the basis values and get some data about the spline spaces // Compute the basis values and get some data about the spline spaces
if (u_from_right) { if (u_from_right) {
surf->basis_u().computeBasisValues(fe.u, &b0[0], derivs, resolution); basis->basis_u().computeBasisValues(fe.u, &b0[0], derivs, resolution);
} else { } else {
surf->basis_u().computeBasisValuesLeft(fe.u, &b0[0], derivs, resolution); basis->basis_u().computeBasisValuesLeft(fe.u, &b0[0], derivs, resolution);
} }
int uleft = surf->basis_u().lastKnotInterval(); int uleft = basis->basis_u().lastKnotInterval();
if (v_from_right) { if (v_from_right) {
surf->basis_v().computeBasisValues(fe.v, &b1[0], derivs, resolution); basis->basis_v().computeBasisValues(fe.v, &b1[0], derivs, resolution);
} else { } else {
surf->basis_v().computeBasisValuesLeft(fe.v, &b1[0], derivs, resolution); basis->basis_v().computeBasisValuesLeft(fe.v, &b1[0], derivs, resolution);
} }
int vleft = surf->basis_v().lastKnotInterval(); int vleft = basis->basis_v().lastKnotInterval();
// Compute the tensor product value // Compute the tensor product value
int coefind = uleft-uorder+1 + unum*(vleft-vorder+1); int coefind = uleft-uorder+1 + unum*(vleft-vorder+1);

14
src/ASM/SplineFields2D.h
View File

@ -33,9 +33,16 @@ class SplineFields2D : public Fields
{ {
public: public:
//! \brief The constructor sets the field name. //! \brief The constructor sets the field name.
//! \param[in] geometry Spline surface geometry //! \param[in] bf Spline basis description
//! \param[in] name Name of spline field //! \param[in] name Name of spline field
SplineFields2D(Go::SplineSurface *geometry, char* name = NULL); 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);
//! \brief Empty destructor. //! \brief Empty destructor.
virtual ~SplineFields2D(); virtual ~SplineFields2D();
@ -68,7 +75,8 @@ public:
bool gradCoor(const Vec3& x, Matrix& grad) const; bool gradCoor(const Vec3& x, Matrix& grad) const;
protected: protected:
Go::SplineSurface* surf; //!< Spline surface geometry description Go::SplineSurface* basis; //!< Spline basis description
Go::SplineSurface* surf; //!< Spline geoemtry description
}; };
#endif #endif

115
src/ASM/SplineFields3D.C
View File

@ -24,17 +24,17 @@ namespace Go {
} }
SplineFields3D::SplineFields3D(Go::SplineVolume *geometry, char* name) SplineFields3D::SplineFields3D(Go::SplineVolume *bf, char* name)
: Fields(3,name), vol(geometry) : Fields(3,name), basis(bf), vol(bf)
{ {
const int n1 = vol->numCoefs(0); const int n1 = basis->numCoefs(0);
const int n2 = vol->numCoefs(1); const int n2 = basis->numCoefs(1);
const int n3 = vol->numCoefs(2); const int n3 = basis->numCoefs(2);
nno = n1*n2*n3; nno = n1*n2*n3;
const int p1 = vol->order(0); const int p1 = basis->order(0);
const int p2 = vol->order(1); const int p2 = basis->order(1);
const int p3 = vol->order(2); const int p3 = basis->order(2);
nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1); nelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1);
// Number of fields set in fill // Number of fields set in fill
@ -42,10 +42,31 @@ SplineFields3D::SplineFields3D(Go::SplineVolume *geometry, char* name)
} }
SplineFields3D::SplineFields3D(Go::SplineVolume *bf,
Go::SplineVolume *geometry,
char* name)
: Fields(3,name), basis(bf), vol(geometry)
{
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);
// Number of fields set in fill
nf = 0;
}
SplineFields3D::~SplineFields3D() SplineFields3D::~SplineFields3D()
{ {
// Set geometry to NULL; should not be deallocated here // Set geometry to NULL; should not be deallocated here
vol = NULL; basis = NULL;
} }
@ -58,24 +79,24 @@ bool SplineFields3D::valueNode(int node, Vector& vals) const
bool SplineFields3D::valueFE(const FiniteElement& fe, Vector& vals) const bool SplineFields3D::valueFE(const FiniteElement& fe, Vector& vals) const
{ {
if (!vol) return false; if (!basis) return false;
Go::Point pt; Go::Point pt;
// vol->pointValue(pt,values,fe.u,fe.v,fe.w); // basis->pointValue(pt,values,fe.u,fe.v,fe.w);
// vals.resize(pt.size()); // vals.resize(pt.size());
// for (int i = 0;i < pt.size();i++) // for (int i = 0;i < pt.size();i++)
// vals[i] = pt[i]; // vals[i] = pt[i];
const int uorder = vol->order(0); const int uorder = basis->order(0);
const int vorder = vol->order(1); const int vorder = basis->order(1);
const int worder = vol->order(2); const int worder = basis->order(2);
const int unum = vol->numCoefs(0); const int unum = basis->numCoefs(0);
const int vnum = vol->numCoefs(1); const int vnum = basis->numCoefs(1);
const int wnum = vol->numCoefs(2); const int wnum = basis->numCoefs(2);
const int dim = vol->dimension(); const int dim = basis->dimension();
const bool rational = vol->rational(); const bool rational = basis->rational();
int kdim = rational ? dim + 1 : dim; int kdim = rational ? dim + 1 : dim;
const int ncomp = values.size()/(unum*vnum*wnum); const int ncomp = values.size()/(unum*vnum*wnum);
@ -95,12 +116,12 @@ bool SplineFields3D::valueFE(const FiniteElement& fe, Vector& vals) const
tempResult.resize(ncomp); tempResult.resize(ncomp);
// compute tbe basis values and get some data about the spline spaces // compute tbe basis values and get some data about the spline spaces
vol->basis(0).computeBasisValues(fe.u, Bu.begin()); basis->basis(0).computeBasisValues(fe.u, Bu.begin());
vol->basis(1).computeBasisValues(fe.v, Bv.begin()); basis->basis(1).computeBasisValues(fe.v, Bv.begin());
vol->basis(2).computeBasisValues(fe.w, Bw.begin()); basis->basis(2).computeBasisValues(fe.w, Bw.begin());
const int uleft = vol->basis(0).lastKnotInterval(); const int uleft = basis->basis(0).lastKnotInterval();
const int vleft = vol->basis(1).lastKnotInterval(); const int vleft = basis->basis(1).lastKnotInterval();
const int wleft = vol->basis(2).lastKnotInterval(); const int wleft = basis->basis(2).lastKnotInterval();
// compute the tensor product value // compute the tensor product value
const int val_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1 + vnum * (wleft - worder + 1))) * ncomp; const int val_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1 + vnum * (wleft - worder + 1))) * ncomp;
@ -115,7 +136,7 @@ bool SplineFields3D::valueFE(const FiniteElement& fe, Vector& vals) const
w = 0.0; w = 0.0;
const int w_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1 + vnum * (wleft - worder + 1))) * kdim + dim; const int w_start_ix = (uleft - uorder + 1 + unum * (vleft - vorder + 1 + vnum * (wleft - worder + 1))) * kdim + dim;
const double *w_ptr = &(vol->rcoefs_begin()[w_start_ix]); const double *w_ptr = &(basis->rcoefs_begin()[w_start_ix]);
for (double* bval_w_ptr = Bw.begin(); bval_w_ptr != Bw.end(); ++bval_w_ptr) { for (double* bval_w_ptr = Bw.begin(); bval_w_ptr != Bw.end(); ++bval_w_ptr) {
const double bval_w = *bval_w_ptr; const double bval_w = *bval_w_ptr;
@ -204,11 +225,11 @@ bool SplineFields3D::valueCoor(const Vec3& x, Vector& vals) const
bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const
{ {
if (!vol) return false; if (!basis) return false;
// // Derivatives of solution in parametric space // // Derivatives of solution in parametric space
// std::vector<Go::Point> pts; // std::vector<Go::Point> pts;
// vol->pointValue(pts,values,fe.u,fe.v,fe.w,1); // basis->pointValue(pts,values,fe.u,fe.v,fe.w,1);
// // Gradient of field wrt parametric coordinates // // Gradient of field wrt parametric coordinates
// Matrix gradXi(nf,3); // Matrix gradXi(nf,3);
@ -218,7 +239,7 @@ bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const
// // Derivatives of coordinate mapping // // Derivatives of coordinate mapping
// std::vector<Go::Point> xpts(3); // std::vector<Go::Point> xpts(3);
// vol->point(xpts,fe.u,fe.v,fe.w,1); // basis->point(xpts,fe.u,fe.v,fe.w,1);
// // Jacobian matrix // // Jacobian matrix
// Matrix Jac(3,3); // Matrix Jac(3,3);
@ -237,13 +258,13 @@ bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const
int derivs = 1; int derivs = 1;
int totpts = (derivs + 1)*(derivs + 2)*(derivs + 3)/6; int totpts = (derivs + 1)*(derivs + 2)*(derivs + 3)/6;
const int unum = vol->numCoefs(0); const int unum = basis->numCoefs(0);
const int vnum = vol->numCoefs(1); const int vnum = basis->numCoefs(1);
const int wnum = vol->numCoefs(2); const int wnum = basis->numCoefs(2);
const int uorder = vol->order(0); const int uorder = basis->order(0);
const int vorder = vol->order(1); const int vorder = basis->order(1);
const int worder = vol->order(2); const int worder = basis->order(2);
const bool u_from_right = true; const bool u_from_right = true;
const bool v_from_right = true; const bool v_from_right = true;
@ -252,8 +273,8 @@ bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const
const double resolution = 1.0e-12; const double resolution = 1.0e-12;
const int ncomp = values.size()/(unum*vnum*wnum); const int ncomp = values.size()/(unum*vnum*wnum);
const int dim = vol->dimension(); const int dim = basis->dimension();
const bool rational = vol->rational(); const bool rational = basis->rational();
for (int i = 0; i < totpts; ++i) { for (int i = 0; i < totpts; ++i) {
if (pts[i].dimension() != ncomp) { if (pts[i].dimension() != ncomp) {
@ -262,7 +283,7 @@ bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const
} }
// Take care of the rational case // Take care of the rational case
const std::vector<double>::iterator co = rational ? vol->rcoefs_begin() : vol->coefs_begin(); const std::vector<double>::iterator co = rational ? basis->rcoefs_begin() : basis->coefs_begin();
int kdim = dim + (rational ? 1 : 0); int kdim = dim + (rational ? 1 : 0);
int ndim = ncomp + (rational ? 1 : 0); int ndim = ncomp + (rational ? 1 : 0);
@ -278,23 +299,23 @@ bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const
// Compute the basis values and get some data about the spline spaces // Compute the basis values and get some data about the spline spaces
if (u_from_right) if (u_from_right)
vol->basis(0).computeBasisValues(fe.u, &b0[0], derivs, resolution); basis->basis(0).computeBasisValues(fe.u, &b0[0], derivs, resolution);
else else
vol->basis(1).computeBasisValuesLeft(fe.u, &b0[0], derivs, resolution); basis->basis(1).computeBasisValuesLeft(fe.u, &b0[0], derivs, resolution);
int uleft = vol->basis(0).lastKnotInterval(); int uleft = basis->basis(0).lastKnotInterval();
if (v_from_right) if (v_from_right)
vol->basis(1).computeBasisValues(fe.v, &b1[0], derivs, resolution); basis->basis(1).computeBasisValues(fe.v, &b1[0], derivs, resolution);
else else
vol->basis(1).computeBasisValuesLeft(fe.v, &b1[0], derivs, resolution); basis->basis(1).computeBasisValuesLeft(fe.v, &b1[0], derivs, resolution);
int vleft = vol->basis(1).lastKnotInterval(); int vleft = basis->basis(1).lastKnotInterval();
if (w_from_right) if (w_from_right)
vol->basis(2).computeBasisValues(fe.w, &b2[0], derivs, resolution); basis->basis(2).computeBasisValues(fe.w, &b2[0], derivs, resolution);
else else
vol->basis(2).computeBasisValuesLeft(fe.w, &b2[0], derivs, resolution); basis->basis(2).computeBasisValuesLeft(fe.w, &b2[0], derivs, resolution);
int wleft = vol->basis(2).lastKnotInterval(); int wleft = basis->basis(2).lastKnotInterval();
// Compute the tensor product value // Compute the tensor product value
int coefind = uleft-uorder+1 + unum*(vleft-vorder+1 + vnum*(wleft-worder+1)); int coefind = uleft-uorder+1 + unum*(vleft-vorder+1 + vnum*(wleft-worder+1));

14
src/ASM/SplineFields3D.h
View File

@ -33,9 +33,16 @@ class SplineFields3D : public Fields
{ {
public: public:
//! \brief The constructor sets the field name. //! \brief The constructor sets the field name.
//! \param[in] geometry Spline volume geometry //! \param[in] bf Spline basis description
//! \param[in] name Name of spline field //! \param[in] name Name of spline field
SplineFields3D(Go::SplineVolume *geometry, char* name = NULL); 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);
//! \brief Empty destructor. //! \brief Empty destructor.
virtual ~SplineFields3D(); virtual ~SplineFields3D();
@ -68,7 +75,8 @@ public:
bool gradCoor(const Vec3& x, Matrix& grad) const; bool gradCoor(const Vec3& x, Matrix& grad) const;
protected: protected:
Go::SplineVolume* vol; //!< Spline volume geometry description Go::SplineVolume* basis; //!< Spline basis description
Go::SplineVolume* vol; //!< Spline geometry description
}; };
#endif #endif