From 24b88f6a253c04a886fe82b13f42d424bfacf055 Mon Sep 17 00:00:00 2001 From: rho Date: Tue, 29 Nov 2011 10:42:40 +0000 Subject: [PATCH] Added support for mixed formulation git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1303 e10b68d5-8a6e-419e-a041-bce267b0401d --- src/ASM/SplineField2D.C | 85 +++++++++++++++++------------ src/ASM/SplineField2D.h | 14 ++++- src/ASM/SplineField3D.C | 109 +++++++++++++++++++++---------------- src/ASM/SplineField3D.h | 13 ++++- src/ASM/SplineFields2D.C | 90 ++++++++++++++++++------------ src/ASM/SplineFields2D.h | 14 ++++- src/ASM/SplineFields3D.C | 115 +++++++++++++++++++++++---------------- src/ASM/SplineFields3D.h | 14 ++++- 8 files changed, 278 insertions(+), 176 deletions(-) diff --git a/src/ASM/SplineField2D.C b/src/ASM/SplineField2D.C index a3caf990..4824e021 100644 --- a/src/ASM/SplineField2D.C +++ b/src/ASM/SplineField2D.C @@ -19,23 +19,38 @@ #include "GoTools/geometry/SplineUtils.h" -SplineField2D::SplineField2D(Go::SplineSurface *geometry, char* name) - : Field(2,name), surf(geometry) +SplineField2D::SplineField2D(Go::SplineSurface *bf, char* name) + : Field(2,name), basis(bf), surf(bf) { - const int n1 = surf->numCoefs_u(); - const int n2 = surf->numCoefs_v(); + const int n1 = basis->numCoefs_u(); + const int n2 = basis->numCoefs_v(); nno = n1*n2; - const int p1 = surf->order_u(); - const int p2 = surf->order_v(); + const int p1 = basis->order_u(); + 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); } SplineField2D::~SplineField2D() { - // Set geometry pointer to NULL, should not be deallocated here - surf = NULL; + // Set basis and geometry pointers to NULL, should not be deallocated here + basis = surf = NULL; } @@ -49,16 +64,16 @@ double SplineField2D::valueNode(int node) const double SplineField2D::valueFE(const FiniteElement& fe) const { // Go::Point pt; -// surf->pointValue(pt,values,fe.u,fe.v); +// basis->pointValue(pt,values,fe.u,fe.v); // return pt[0]; - const int uorder = surf->order_u(); - const int vorder = surf->order_v(); - const int unum = surf->numCoefs_u(); - //const int vnum = surf->numCoefs_v(); + const int uorder = basis->order_u(); + const int vorder = basis->order_v(); + const int unum = basis->numCoefs_u(); + //const int vnum = basis->numCoefs_v(); - const int dim = surf->dimension(); - const bool rational = surf->rational(); + const int dim = basis->dimension(); + const bool rational = basis->rational(); const int kdim = rational ? dim + 1 : dim; //const int ncomp = values.size()/(unum*vnum); @@ -71,11 +86,11 @@ double SplineField2D::valueFE(const FiniteElement& fe) const Bu.resize(uorder); Bv.resize(vorder); - surf->basis_u().computeBasisValues(fe.u,Bu.begin()); - surf->basis_v().computeBasisValues(fe.v,Bv.begin()); + basis->basis_u().computeBasisValues(fe.u,Bu.begin()); + basis->basis_v().computeBasisValues(fe.v,Bv.begin()); - const int uleft = surf->basis_u().lastKnotInterval(); - const int vleft = surf->basis_v().lastKnotInterval(); + const int uleft = basis->basis_u().lastKnotInterval(); + const int vleft = basis->basis_v().lastKnotInterval(); // compute the tensor product value 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 w = 0.0; 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) { 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 { - if (!surf) return false; + if (!basis) return false; // // Derivatives of solution in parametric space // std::vector 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 // Vector gradXi(2); @@ -149,7 +164,7 @@ bool SplineField2D::gradFE(const FiniteElement& fe, Vector& grad) const // // Derivatives of coordinate mapping // std::vector xpts(3); -// surf->point(xpts,fe.u,fe.v,1); +// basis->point(xpts,fe.u,fe.v,1); // // Jacobian matrix // Matrix Jac(2,2); @@ -167,9 +182,9 @@ bool SplineField2D::gradFE(const FiniteElement& fe, Vector& grad) const // Gradient of field wrt parametric coordinates Vector gradXi(2); - const int uorder = surf->order_u(); - const int vorder = surf->order_v(); - const int unum = surf->numCoefs_u(); + const int uorder = basis->order_u(); + const int vorder = basis->order_v(); + const int unum = basis->numCoefs_u(); const int derivs = 1; const int totpts = 3; @@ -178,11 +193,11 @@ bool SplineField2D::gradFE(const FiniteElement& fe, Vector& grad) const const double resolution = 1.0e-12; // Dimension - const int dim = surf->dimension(); - const bool rational = surf->rational(); + const int dim = basis->dimension(); + const bool rational = basis->rational(); // Take care of the rational case - const std::vector::iterator co = rational ? surf->rcoefs_begin() : surf->coefs_begin(); + const std::vector::iterator co = rational ? basis->rcoefs_begin() : basis->coefs_begin(); int kdim = dim + (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 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 - surf->basis_u().computeBasisValuesLeft(fe.u, &b0[0], derivs, resolution); - int uleft = surf->basis_u().lastKnotInterval(); + basis->basis_u().computeBasisValuesLeft(fe.u, &b0[0], derivs, resolution); + int uleft = basis->basis_u().lastKnotInterval(); 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 - surf->basis_v().computeBasisValuesLeft(fe.v, &b1[0], derivs, resolution); - int vleft = surf->basis_v().lastKnotInterval(); + basis->basis_v().computeBasisValuesLeft(fe.v, &b1[0], derivs, resolution); + int vleft = basis->basis_v().lastKnotInterval(); // Compute the tensor product value int coefind = uleft-uorder+1 + unum*(vleft-vorder+1); diff --git a/src/ASM/SplineField2D.h b/src/ASM/SplineField2D.h index d7074ba1..03ebe0f1 100644 --- a/src/ASM/SplineField2D.h +++ b/src/ASM/SplineField2D.h @@ -33,9 +33,16 @@ class SplineField2D : public Field { public: //! \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 - 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. virtual ~SplineField2D(); @@ -65,7 +72,8 @@ public: bool gradCoor(const Vec3& x, Vector& grad) const; protected: - Go::SplineSurface* surf; //!< Spline surface geometry description + Go::SplineSurface* basis; //!< Spline basis description + Go::SplineSurface* surf; //!< Spline geometry description }; #endif diff --git a/src/ASM/SplineField3D.C b/src/ASM/SplineField3D.C index 2130fdaa..8ea9d665 100644 --- a/src/ASM/SplineField3D.C +++ b/src/ASM/SplineField3D.C @@ -23,25 +23,42 @@ namespace Go { } -SplineField3D::SplineField3D(Go::SplineVolume *geometry, char* name) - : Field(3,name), vol(geometry) +SplineField3D::SplineField3D(Go::SplineVolume *bf, char* name) + : Field(3,name), basis(bf), vol(bf) { - const int n1 = vol->numCoefs(0); - const int n2 = vol->numCoefs(1); - const int n3 = vol->numCoefs(2); + 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 = vol->order(0); - const int p2 = vol->order(1); - const int p3 = vol->order(2); + 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); +} + + +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); } SplineField3D::~SplineField3D() { - // Set geometry pointer to NULL, should not be deallocated here - vol = NULL; + // Set basis and geometry pointers to NULL, should not be deallocated here + basis = vol = NULL; } @@ -55,19 +72,19 @@ double SplineField3D::valueNode(int node) const double SplineField3D::valueFE(const FiniteElement& fe) const { // 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]; double val = 0.0; - const int uorder = vol->order(0); - const int vorder = vol->order(1); - const int worder = vol->order(2); - const int unum = vol->numCoefs(0); - const int vnum = vol->numCoefs(1); + const int uorder = basis->order(0); + const int vorder = basis->order(1); + const int worder = basis->order(2); + const int unum = basis->numCoefs(0); + const int vnum = basis->numCoefs(1); - const int dim = vol->dimension(); - const bool rational = vol->rational(); + const int dim = basis->dimension(); + const bool rational = basis->rational(); int kdim = rational ? dim + 1 : dim; static Go::ScratchVect Bu(uorder); @@ -80,12 +97,12 @@ double SplineField3D::valueFE(const FiniteElement& fe) const Bw.resize(worder); // compute tbe basis values and get some data about the spline spaces - vol->basis(0).computeBasisValues(fe.u, Bu.begin()); - vol->basis(1).computeBasisValues(fe.v, Bv.begin()); - vol->basis(2).computeBasisValues(fe.w, Bw.begin()); - const int uleft = vol->basis(0).lastKnotInterval(); - const int vleft = vol->basis(1).lastKnotInterval(); - const int wleft = vol->basis(2).lastKnotInterval(); + basis->basis(0).computeBasisValues(fe.u, Bu.begin()); + basis->basis(1).computeBasisValues(fe.v, Bv.begin()); + basis->basis(2).computeBasisValues(fe.w, Bw.begin()); + const int uleft = basis->basis(0).lastKnotInterval(); + const int vleft = basis->basis(1).lastKnotInterval(); + const int wleft = basis->basis(2).lastKnotInterval(); // compute the tensor product value 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; 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) { 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 { - if (!vol) return false; + if (!basis) return false; // // Derivatives of solution in parametric space // std::vector 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 // Vector gradXi(3); @@ -176,7 +193,7 @@ bool SplineField3D::gradFE(const FiniteElement& fe, Vector& grad) const // // Derivatives of coordinate mapping // std::vector xpts(3); -// vol->point(xpts,fe.u,fe.v,fe.w,1); +// basis->point(xpts,fe.u,fe.v,fe.w,1); // // Jacobian matrix // Matrix Jac(3,3); @@ -197,13 +214,13 @@ bool SplineField3D::gradFE(const FiniteElement& fe, Vector& grad) const int derivs = 1; int totpts = (derivs + 1)*(derivs + 2)*(derivs + 3)/6; - const int unum = vol->numCoefs(0); - const int vnum = vol->numCoefs(1); - const int wnum = vol->numCoefs(2); + const int unum = basis->numCoefs(0); + const int vnum = basis->numCoefs(1); + const int wnum = basis->numCoefs(2); - const int uorder = vol->order(0); - const int vorder = vol->order(1); - const int worder = vol->order(2); + const int uorder = basis->order(0); + const int vorder = basis->order(1); + const int worder = basis->order(2); const bool u_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 int ncomp = values.size()/(unum*vnum*wnum); - const int dim = vol->dimension(); - const bool rational = vol->rational(); + const int dim = basis->dimension(); + const bool rational = basis->rational(); // Take care of the rational case - const std::vector::iterator co = rational ? vol->rcoefs_begin() : vol->coefs_begin(); + const std::vector::iterator co = rational ? basis->rcoefs_begin() : basis->coefs_begin(); int kdim = dim + (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 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 - 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) - vol->basis(1).computeBasisValues(fe.v, &b1[0], derivs, resolution); + basis->basis(1).computeBasisValues(fe.v, &b1[0], derivs, resolution); 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) - vol->basis(2).computeBasisValues(fe.w, &b2[0], derivs, resolution); + basis->basis(2).computeBasisValues(fe.w, &b2[0], derivs, resolution); 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 int coefind = uleft-uorder+1 + unum*(vleft-vorder+1 + vnum*(wleft-worder+1)); diff --git a/src/ASM/SplineField3D.h b/src/ASM/SplineField3D.h index 53a72489..c4bd570a 100644 --- a/src/ASM/SplineField3D.h +++ b/src/ASM/SplineField3D.h @@ -33,9 +33,15 @@ class SplineField3D : public Field { public: //! \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 - 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. virtual ~SplineField3D(); @@ -65,7 +71,8 @@ public: bool gradCoor(const Vec3& x, Vector& grad) const; protected: - Go::SplineVolume* vol; //!< Spline volume geometry description + Go::SplineVolume* basis; //!< Spline basis description + Go::SplineVolume* vol; //!< Spline volume description }; #endif diff --git a/src/ASM/SplineFields2D.C b/src/ASM/SplineFields2D.C index 42438355..23d4e8c0 100644 --- a/src/ASM/SplineFields2D.C +++ b/src/ASM/SplineFields2D.C @@ -22,15 +22,33 @@ #endif -SplineFields2D::SplineFields2D(Go::SplineSurface *geometry, char* name) - : Fields(2,name), surf(geometry) +SplineFields2D::SplineFields2D(Go::SplineSurface *bf, char* name) + : Fields(2,name), basis(bf), surf(bf) { - const int n1 = surf->numCoefs_u(); - const int n2 = surf->numCoefs_v(); + const int n1 = basis->numCoefs_u(); + const int n2 = basis->numCoefs_v(); nno = n1*n2; - const int p1 = surf->order_u(); - const int p2 = surf->order_v(); + 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; +} + + +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); // Number of fields set in fill @@ -41,7 +59,7 @@ SplineFields2D::SplineFields2D(Go::SplineSurface *geometry, char* name) SplineFields2D::~SplineFields2D() { // 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 { - if (!surf) return false; + if (!basis) return false; // Go::Point pt; -// surf->pointValue(pt,values,fe.u,fe.v); +// basis->pointValue(pt,values,fe.u,fe.v); // vals.resize(pt.size()); // for (int i = 0;i < pt.size();i++) // vals[i] = pt[i]; - const int uorder = surf->order_u(); - const int vorder = surf->order_v(); - const int unum = surf->numCoefs_u(); - const int vnum = surf->numCoefs_v(); + const int uorder = basis->order_u(); + const int vorder = basis->order_v(); + const int unum = basis->numCoefs_u(); + const int vnum = basis->numCoefs_v(); - const int dim = surf->dimension(); - const bool rational = surf->rational(); + const int dim = basis->dimension(); + const bool rational = basis->rational(); int kdim = rational ? dim + 1 : dim; const int ncomp = values.size()/(unum*vnum); @@ -83,11 +101,11 @@ bool SplineFields2D::valueFE(const FiniteElement& fe, Vector& vals) const Bv.resize(vorder); tempVal.resize(ncomp); - surf->basis_u().computeBasisValues(fe.u,Bu.begin()); - surf->basis_v().computeBasisValues(fe.v,Bv.begin()); + basis->basis_u().computeBasisValues(fe.u,Bu.begin()); + basis->basis_v().computeBasisValues(fe.v,Bv.begin()); - const int uleft = surf->basis_u().lastKnotInterval(); - const int vleft = surf->basis_v().lastKnotInterval(); + const int uleft = basis->basis_u().lastKnotInterval(); + const int vleft = basis->basis_v().lastKnotInterval(); //compute the tensor product value 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; 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) { 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 { - if (!surf) return false; + if (!basis) return false; // // Derivatives of solution in parametric space // std::vector 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 // Matrix gradXi(nf,2); @@ -183,7 +201,7 @@ bool SplineFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const // // Derivatives of coordinate mapping // std::vector xpts(3); -// surf->point(xpts,fe.u,fe.v,1); +// basis->point(xpts,fe.u,fe.v,1); // // Jacobian matrix // Matrix Jac(2,2); @@ -201,10 +219,10 @@ bool SplineFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const // Gradient of field Matrix gradXi(nf,2); - const int uorder = surf->order_u(); - const int vorder = surf->order_v(); - const int unum = surf->numCoefs_u(); - const int vnum = surf->numCoefs_v(); + const int uorder = basis->order_u(); + const int vorder = basis->order_v(); + const int unum = basis->numCoefs_u(); + const int vnum = basis->numCoefs_v(); const int derivs = 1; 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); // Dimension - const int dim = surf->dimension(); - const bool rational = surf->rational(); + const int dim = basis->dimension(); + const bool rational = basis->rational(); // Take care of the rational case - const std::vector::iterator co = rational ? surf->rcoefs_begin() : surf->coefs_begin(); + const std::vector::iterator co = rational ? basis->rcoefs_begin() : basis->coefs_begin(); int kdim = dim + (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 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 { - 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) { - surf->basis_v().computeBasisValues(fe.v, &b1[0], derivs, resolution); + basis->basis_v().computeBasisValues(fe.v, &b1[0], derivs, resolution); } 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 int coefind = uleft-uorder+1 + unum*(vleft-vorder+1); diff --git a/src/ASM/SplineFields2D.h b/src/ASM/SplineFields2D.h index df0bb103..83b1d4cf 100644 --- a/src/ASM/SplineFields2D.h +++ b/src/ASM/SplineFields2D.h @@ -33,9 +33,16 @@ class SplineFields2D : public Fields { public: //! \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 - 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. virtual ~SplineFields2D(); @@ -68,7 +75,8 @@ public: bool gradCoor(const Vec3& x, Matrix& grad) const; protected: - Go::SplineSurface* surf; //!< Spline surface geometry description + Go::SplineSurface* basis; //!< Spline basis description + Go::SplineSurface* surf; //!< Spline geoemtry description }; #endif diff --git a/src/ASM/SplineFields3D.C b/src/ASM/SplineFields3D.C index 32dc2e46..dc039edb 100644 --- a/src/ASM/SplineFields3D.C +++ b/src/ASM/SplineFields3D.C @@ -24,17 +24,17 @@ namespace Go { } -SplineFields3D::SplineFields3D(Go::SplineVolume *geometry, char* name) - : Fields(3,name), vol(geometry) +SplineFields3D::SplineFields3D(Go::SplineVolume *bf, char* name) + : Fields(3,name), basis(bf), vol(bf) { - const int n1 = vol->numCoefs(0); - const int n2 = vol->numCoefs(1); - const int n3 = vol->numCoefs(2); + 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 = vol->order(0); - const int p2 = vol->order(1); - const int p3 = vol->order(2); + 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 @@ -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() { // 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 { - if (!vol) return false; + if (!basis) return false; 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()); // for (int i = 0;i < pt.size();i++) // vals[i] = pt[i]; - const int uorder = vol->order(0); - const int vorder = vol->order(1); - const int worder = vol->order(2); - const int unum = vol->numCoefs(0); - const int vnum = vol->numCoefs(1); - const int wnum = vol->numCoefs(2); + const int uorder = basis->order(0); + const int vorder = basis->order(1); + const int worder = basis->order(2); + const int unum = basis->numCoefs(0); + const int vnum = basis->numCoefs(1); + const int wnum = basis->numCoefs(2); - const int dim = vol->dimension(); - const bool rational = vol->rational(); + const int dim = basis->dimension(); + const bool rational = basis->rational(); int kdim = rational ? dim + 1 : dim; const int ncomp = values.size()/(unum*vnum*wnum); @@ -95,12 +116,12 @@ bool SplineFields3D::valueFE(const FiniteElement& fe, Vector& vals) const tempResult.resize(ncomp); // compute tbe basis values and get some data about the spline spaces - vol->basis(0).computeBasisValues(fe.u, Bu.begin()); - vol->basis(1).computeBasisValues(fe.v, Bv.begin()); - vol->basis(2).computeBasisValues(fe.w, Bw.begin()); - const int uleft = vol->basis(0).lastKnotInterval(); - const int vleft = vol->basis(1).lastKnotInterval(); - const int wleft = vol->basis(2).lastKnotInterval(); + basis->basis(0).computeBasisValues(fe.u, Bu.begin()); + basis->basis(1).computeBasisValues(fe.v, Bv.begin()); + basis->basis(2).computeBasisValues(fe.w, Bw.begin()); + const int uleft = basis->basis(0).lastKnotInterval(); + const int vleft = basis->basis(1).lastKnotInterval(); + const int wleft = basis->basis(2).lastKnotInterval(); // compute the tensor product value 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; 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) { 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 { - if (!vol) return false; + if (!basis) return false; // // Derivatives of solution in parametric space // std::vector 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 // Matrix gradXi(nf,3); @@ -218,7 +239,7 @@ bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const // // Derivatives of coordinate mapping // std::vector xpts(3); -// vol->point(xpts,fe.u,fe.v,fe.w,1); +// basis->point(xpts,fe.u,fe.v,fe.w,1); // // Jacobian matrix // Matrix Jac(3,3); @@ -237,13 +258,13 @@ bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const int derivs = 1; int totpts = (derivs + 1)*(derivs + 2)*(derivs + 3)/6; - const int unum = vol->numCoefs(0); - const int vnum = vol->numCoefs(1); - const int wnum = vol->numCoefs(2); + const int unum = basis->numCoefs(0); + const int vnum = basis->numCoefs(1); + const int wnum = basis->numCoefs(2); - const int uorder = vol->order(0); - const int vorder = vol->order(1); - const int worder = vol->order(2); + const int uorder = basis->order(0); + const int vorder = basis->order(1); + const int worder = basis->order(2); const bool u_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 int ncomp = values.size()/(unum*vnum*wnum); - const int dim = vol->dimension(); - const bool rational = vol->rational(); + const int dim = basis->dimension(); + const bool rational = basis->rational(); for (int i = 0; i < totpts; ++i) { if (pts[i].dimension() != ncomp) { @@ -262,7 +283,7 @@ bool SplineFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const } // Take care of the rational case - const std::vector::iterator co = rational ? vol->rcoefs_begin() : vol->coefs_begin(); + const std::vector::iterator co = rational ? basis->rcoefs_begin() : basis->coefs_begin(); int kdim = dim + (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 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 - 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) - vol->basis(1).computeBasisValues(fe.v, &b1[0], derivs, resolution); + basis->basis(1).computeBasisValues(fe.v, &b1[0], derivs, resolution); 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) - vol->basis(2).computeBasisValues(fe.w, &b2[0], derivs, resolution); + basis->basis(2).computeBasisValues(fe.w, &b2[0], derivs, resolution); 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 int coefind = uleft-uorder+1 + unum*(vleft-vorder+1 + vnum*(wleft-worder+1)); diff --git a/src/ASM/SplineFields3D.h b/src/ASM/SplineFields3D.h index 28479a8d..e8466659 100644 --- a/src/ASM/SplineFields3D.h +++ b/src/ASM/SplineFields3D.h @@ -33,9 +33,16 @@ class SplineFields3D : public Fields { public: //! \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 - 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. virtual ~SplineFields3D(); @@ -68,7 +75,8 @@ public: bool gradCoor(const Vec3& x, Matrix& grad) const; protected: - Go::SplineVolume* vol; //!< Spline volume geometry description + Go::SplineVolume* basis; //!< Spline basis description + Go::SplineVolume* vol; //!< Spline geometry description }; #endif