From 58e55f5b88165b484795bbdc3ed7a9c4aa46ed63 Mon Sep 17 00:00:00 2001 From: rho Date: Tue, 10 Jan 2012 12:09:56 +0000 Subject: [PATCH] Made getElementCoordinates and getNodalCoordinates public. Added function getNoElmDOF. Added some functionality for computation of stabilization parameters git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1363 e10b68d5-8a6e-419e-a041-bce267b0401d --- src/ASM/ASMbase.h | 17 +++++++++++++---- src/ASM/ASMmxBase.C | 17 +++++++++++++++-- src/ASM/ASMmxBase.h | 11 +++++++++-- src/ASM/ASMs1D.h | 21 +++++++++++---------- src/ASM/ASMs1DLag.h | 21 +++++++++++---------- src/ASM/ASMs2D.C | 42 +++++++++++++++++++++++++++++++++++------- src/ASM/ASMs2D.h | 20 +++++++++++--------- src/ASM/ASMs2DLag.h | 28 +++++++++++----------------- src/ASM/ASMs2Dmx.C | 21 ++++++++++++++++----- src/ASM/ASMs2Dmx.h | 12 ++++++------ src/ASM/ASMs3D.C | 33 ++++++++++++++++++++++++++++++++- src/ASM/ASMs3D.h | 21 +++++++++++---------- src/ASM/ASMs3DLag.h | 6 ------ src/ASM/ASMs3Dmx.C | 16 ++++++++++++++++ src/ASM/ASMs3Dmx.h | 12 ++++++------ 15 files changed, 203 insertions(+), 95 deletions(-) diff --git a/src/ASM/ASMbase.h b/src/ASM/ASMbase.h index e28bf908..50729fb3 100644 --- a/src/ASM/ASMbase.h +++ b/src/ASM/ASMbase.h @@ -114,6 +114,8 @@ public: unsigned char getNoParamDim() const { return ndim; } //! \brief Returns the number of solution fields. virtual unsigned char getNoFields(int b = 0) const { return b > 1 ? 0 : nf; } + //! \brief Returns the number of solution fields. + virtual size_t getNoElmDOF(int b = 0) const { return b > 1 ? 0 : neldof; } //! \brief Returns local 1-based index of the node with given global number. //! \details If the given node number is not present, 0 is returned. @@ -136,6 +138,12 @@ public: //! \param[out] X nsd\f$\times\f$n-matrix, where \a n is the number of nodes //! in the patch virtual void getNodalCoordinates(Matrix& X) const = 0; + //! \brief Returns a matrix with nodal coordinates for an element. + //! \param[in] iel Element index + //! \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 = 0; + //! \brief Prints out the nodal coordinates of this patch to the given stream. void printNodes(std::ostream& os, const char* heading = 0) const; @@ -440,10 +448,11 @@ public: protected: // Standard finite element data structures - unsigned char ndim; //!< Number of parametric dimensions (1, 2 or 3) - unsigned char nsd; //!< Number of space dimensions (ndim <= nsd <= 3) - unsigned char nf; //!< Number of primary solution fields (1 or larger) - + unsigned char ndim; //!< Number of parametric dimensions (1, 2 or 3) + unsigned char nsd; //!< Number of space dimensions (ndim <= nsd <= 3) + unsigned char nf; //!< Number of primary solution fields (1 or larger) + size_t neldof; //!< Number of degrees of freedom per element + const IntVec& MLGE; //!< Matrix of Local to Global Element numbers const IntVec& MLGN; //!< Matrix of Local to Global Node numbers const IntMat& MNPC; //!< Matrix of Nodal Point Correspondance diff --git a/src/ASM/ASMmxBase.C b/src/ASM/ASMmxBase.C index 13b60fc3..2127f1e0 100644 --- a/src/ASM/ASMmxBase.C +++ b/src/ASM/ASMmxBase.C @@ -14,8 +14,9 @@ #include "ASMmxBase.h" -bool ASMmxBase::geoUsesBasis1 = false; -bool ASMmxBase::useCpminus1 = false; +bool ASMmxBase::geoUsesBasis1 = false; +bool ASMmxBase::useCpminus1 = false; +bool ASMmxBase::useLowOrderBasis1 = false; ASMmxBase::ASMmxBase (unsigned char n_f1, unsigned char n_f2) @@ -25,6 +26,18 @@ ASMmxBase::ASMmxBase (unsigned char n_f1, unsigned char n_f2) } +size_t ASMmxBase::getNoElmDOF(int basis) const +{ + switch(basis) { + case 1: return neldof1; + case 2: return neldof2; + } + + + return neldof1 + neldof2; +} + + void ASMmxBase::initMx (const std::vector& MLGN, const int* sysMadof) { MADOF.resize(MLGN.size()); diff --git a/src/ASM/ASMmxBase.h b/src/ASM/ASMmxBase.h index 4e8b3b4c..9ffee93f 100644 --- a/src/ASM/ASMmxBase.h +++ b/src/ASM/ASMmxBase.h @@ -49,8 +49,12 @@ protected: const std::vector& nodes) const; public: - static bool geoUsesBasis1; //!< If \e true, 1st basis represents the geometry - static bool useCpminus1; //!< If \e true, enforce \f$C^{p-1}\f$ continuity + static bool geoUsesBasis1; //!< If \e true, 1st basis represents the geometry + static bool useCpminus1; //!< If \e true, enforce \f$C^{p-1}\f$ continuity + static bool useLowOrderBasis1; //!< If \e true, basis 1 is of lowest order + + //! \brief Returns the number of solution fields. + virtual size_t getNoElmDOF(int b = 0) const; private: std::vector MADOF; //!< Matrix of accumulated DOFs for this patch @@ -61,6 +65,9 @@ protected: unsigned char nf1; //!< Number of solution fields using first basis unsigned char nf2; //!< Number of solution fields using second basis + + size_t neldof1; //!< Number of degrees of freedom per element for basis 1 + size_t neldof2; //!< Number of degrees of freedom per element for basis 2 }; #endif diff --git a/src/ASM/ASMs1D.h b/src/ASM/ASMs1D.h index 56cd06af..6d68a629 100644 --- a/src/ASM/ASMs1D.h +++ b/src/ASM/ASMs1D.h @@ -59,6 +59,17 @@ public: //! \param[in] inod 1-based node index local to current patch virtual Vec3 getCoord(size_t inod) const; + //! \brief Returns a matrix with nodal coordinates for an element. + //! \param[in] iel Element index + //! \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; + + //! \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 + virtual void getNodalCoordinates(Matrix& X) const; + //! \brief Updates the nodal coordinates for this patch. //! \param[in] displ Incremental displacements to update the coordinates with virtual bool updateCoords(const Vector& displ); @@ -223,16 +234,6 @@ protected: //! \param[in] iel Element index double getParametricLength(int iel) const; - //! \brief Returns a matrix with nodal coordinates for an element. - //! \param[in] iel Element index - //! \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; - //! \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 - virtual void getNodalCoordinates(Matrix& X) const; - //! \brief Returns the number of nodal points in the patch. //! \param[in] basis Which basis to return size parameters for (mixed methods) virtual int getSize(int basis = 0) const; diff --git a/src/ASM/ASMs1DLag.h b/src/ASM/ASMs1DLag.h index 2832cfb8..01792b5c 100644 --- a/src/ASM/ASMs1DLag.h +++ b/src/ASM/ASMs1DLag.h @@ -47,6 +47,17 @@ public: //! This is used to reinitialize the patch after it has been refined. virtual void clear(bool retainGeometry = false); + //! \brief Returns a matrix with nodal coordinates for an element. + //! \param[in] iel Element index + //! \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; + + //! \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 + virtual void getNodalCoordinates(Matrix& X) const; + //! \brief Returns the global coordinates for the given node. //! \param[in] inod 1-based node index local to current patch virtual Vec3 getCoord(size_t inod) const; @@ -131,16 +142,6 @@ protected: // Internal utility methods // ======================== - //! \brief Returns a matrix with nodal coordinates for an element. - //! \param[in] iel Element index - //! \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; - //! \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 - virtual void getNodalCoordinates(Matrix& X) const; - //! \brief Returns the number of nodal points in the patch. virtual int getSize(int = 0) const { return nx; } diff --git a/src/ASM/ASMs2D.C b/src/ASM/ASMs2D.C index 428058e0..24fd4519 100644 --- a/src/ASM/ASMs2D.C +++ b/src/ASM/ASMs2D.C @@ -86,6 +86,9 @@ bool ASMs2D::read (std::istream& is) nsd = surf->dimension(); } + // Define number of DOFs pr element + neldof = surf->order_u()*surf->order_v()*nf; + geo = surf; return true; } @@ -668,6 +671,7 @@ double ASMs2D::getParametricArea (int iel) const double du = surf->knotSpan(0,nodeInd[inod1].I); double dv = surf->knotSpan(1,nodeInd[inod1].J); + return du*dv; } @@ -889,6 +893,25 @@ static double getElmSize (int p1, int p2, const Matrix& X) } +/*! + \brief Computes the characteristic element length from inverse Jacobian. +*/ + +static bool getG (const Matrix& Jinv, const Vector& du, Matrix& G) +{ + const size_t nsd = Jinv.cols(); + + G.resize(nsd,nsd,true); + + for (size_t k = 1;k <= nsd;k++) + for (size_t l = 1;l <= nsd;l++) + for (size_t m = 1;m <= nsd;m++) + G(k,l) += Jinv(m,k)*Jinv(m,l)/(du(k)*du(l)); + + return true; +} + + void ASMs2D::extractBasis (const Go::BasisDerivsSf& spline, Vector& N, Matrix& dNdu) { @@ -994,7 +1017,6 @@ bool ASMs2D::integrate (Integrand& integrand, Matrix3D d2Ndu2, Hess; Vec4 X; - // === Assembly loop over all elements in the patch ========================== int iel = 1; @@ -1012,9 +1034,8 @@ bool ASMs2D::integrate (Integrand& integrand, if (!this->getElementCoordinates(Xnod,iel)) return false; // Compute characteristic element length, if needed - if (integrand.getIntegrandType() == 2) - fe.h = getElmSize(p1,p2,Xnod); - + if (integrand.getIntegrandType() == 2) + fe.h = getElmSize(p1,p2,Xnod); else if (integrand.getIntegrandType() == 3) { // --- Compute average value of basis functions over the element ------- @@ -1062,7 +1083,6 @@ bool ASMs2D::integrate (Integrand& integrand, // If the array is shorter than this, expect a segmentation fault. LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1]; - if (integrand.getIntegrandType() > 10) { // --- Selective reduced integration loop ------------------------------ @@ -1121,10 +1141,18 @@ bool ASMs2D::integrate (Integrand& integrand, if (fe.detJxW == 0.0) continue; // skip singular points // Compute Hessian of coordinate mapping and 2nd order derivatives - if (integrand.getIntegrandType() == 2) + if (integrand.getIntegrandType() == 2) { if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu)) return false; + // RUNAR + Vector dXidu(2); + dXidu(1) = this->getParametricLength(iel,1); + dXidu(2) = this->getParametricLength(iel,2); + if (!getG(Jac,dXidu,fe.G)) + return false; + } + #if SP_DEBUG > 4 std::cout <<"\niel, ip = "<< iel <<" "<< ip <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX << std::endl; @@ -1146,7 +1174,7 @@ bool ASMs2D::integrate (Integrand& integrand, // Assembly of global system integral if (!glInt.assemble(elmInt,fe.iel)) - return false; + return false; } return true; diff --git a/src/ASM/ASMs2D.h b/src/ASM/ASMs2D.h index 38beec94..03f282bc 100644 --- a/src/ASM/ASMs2D.h +++ b/src/ASM/ASMs2D.h @@ -100,6 +100,17 @@ public: //! This is used to reinitialize the patch after it has been refined. virtual void clear(bool retainGeometry = false); + //! \brief Returns a matrix with nodal coordinates for an element. + //! \param[in] iel Element index + //! \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; + + //! \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 + virtual void getNodalCoordinates(Matrix& X) const; + //! \brief Returns the global coordinates for the given node. //! \param[in] inod 1-based node index local to current patch virtual Vec3 getCoord(size_t inod) const; @@ -331,15 +342,6 @@ protected: //! \param[in] iel Element index //! \param[in] dir Local index of the boundary edge double getParametricLength(int iel, int dir) const; - //! \brief Returns a matrix with nodal coordinates for an element. - //! \param[in] iel Element index - //! \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; - //! \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 - virtual void getNodalCoordinates(Matrix& X) const; //! \brief Returns the number of nodal points in each parameter direction. //! \param[out] n1 Number of nodes in first (u) direction diff --git a/src/ASM/ASMs2DLag.h b/src/ASM/ASMs2DLag.h index c49ada67..03c216dd 100644 --- a/src/ASM/ASMs2DLag.h +++ b/src/ASM/ASMs2DLag.h @@ -47,6 +47,17 @@ public: //! This is used to reinitialize the patch after it has been refined. virtual void clear(bool retainGeometry = false); + //! \brief Returns a matrix with nodal coordinates for an element. + //! \param[in] iel Element index + //! \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; + + //! \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 + virtual void getNodalCoordinates(Matrix& X) const; + //! \brief Returns the global coordinates for the given node. //! \param[in] inod 1-based node index local to current patch virtual Vec3 getCoord(size_t inod) const; @@ -130,23 +141,6 @@ public: virtual bool evalSolution(Matrix& sField, const Integrand& integrand, const RealArray* gpar, bool regular = true) const; -protected: - - // Internal utility methods - // ======================== - - //! \brief Returns a matrix with nodal coordinates for an element. - //! \param[in] iel Element index - //! \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 - virtual void getNodalCoordinates(Matrix& X) const; - //! \brief Returns the number of nodal points in each parameter direction. //! \param[out] n1 Number of nodes in first (u) direction //! \param[out] n2 Number of nodes in second (v) direction diff --git a/src/ASM/ASMs2Dmx.C b/src/ASM/ASMs2Dmx.C index 54fc953f..082f4ad5 100644 --- a/src/ASM/ASMs2Dmx.C +++ b/src/ASM/ASMs2Dmx.C @@ -189,7 +189,7 @@ bool ASMs2Dmx::generateFEMTopology () ug,vg,XYZ,ndim, false,XYZ); } - else + else { // Order-elevate basis1 such that it is of one degree higher than basis2 // but only C^p-2 continuous @@ -197,6 +197,12 @@ bool ASMs2Dmx::generateFEMTopology () basis1->raiseOrder(1,1); } basis2 = surf; + + if (useLowOrderBasis1) { + basis2 = basis1; + basis1 = surf; + surf = basis2; + } // Define which basis that should be used to represent the geometry if (geoUsesBasis1) surf = basis1; @@ -249,7 +255,7 @@ bool ASMs2Dmx::generateFEMTopology () if (p1 > n1 || p2 > n2) return false; if (q1 > m1 || q2 > m2) return false; - myMLGE.resize((n1-p1+1)*(n2-p2+1),0); + myMLGE.resize(geoUsesBasis1 ? (n1-p1+1)*(n2-p2+1) : (m1-q1+1)*(m2-q2+1),0); myMLGN.resize(nb1 + nb2); myMNPC.resize(myMLGE.size()); myNodeInd.resize(myMLGN.size()); @@ -340,7 +346,7 @@ bool ASMs2Dmx::generateFEMTopology () iel = inod = 0; for (i2 = 1; i2 <= n2; i2++) for (i1 = 1; i1 <= n1; i1++, inod++) - if (i1 >= p1 && i2 >= p2) + if (i1 >= p1 && i2 >= p2) { if (basis1->knotSpan(0,i1-1) > 0.0) if (basis1->knotSpan(1,i2-1) > 0.0) { @@ -352,9 +358,15 @@ bool ASMs2Dmx::generateFEMTopology () myMNPC[iel][lnod++] = inod - n1*j2 - j1; iel++; - } + } + } } + // Number of DOFs pr element for each basis + neldof1 = basis1->order_u()*basis1->order_v()*nf1; + neldof2 = basis2->order_u()*basis2->order_v()*nf2; + neldof = neldof1 + neldof2; + #ifdef SP_DEBUG std::cout <<"NEL = "<< MLGE.size() <<" NNOD = "<< MLGN.size() << std::endl; #endif @@ -499,7 +511,6 @@ bool ASMs2Dmx::integrate (Integrand& integrand, Matrix dN1du, dN2du, Xnod, Jac; Vec4 X; - // === Assembly loop over all elements in the patch ========================== int iel = 1; diff --git a/src/ASM/ASMs2Dmx.h b/src/ASM/ASMs2Dmx.h index 35e6e094..4dc9b7b1 100644 --- a/src/ASM/ASMs2Dmx.h +++ b/src/ASM/ASMs2Dmx.h @@ -60,6 +60,12 @@ public: //! This is used to reinitialize the patch after it has been refined. virtual void clear(bool retainGeometry = false); + //! \brief Returns a matrix with nodal coordinates for an element. + //! \param[in] iel Element index + //! \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; + //! \brief Returns the global coordinates for the given node. //! \param[in] inod 1-based node index local to current patch virtual Vec3 getCoord(size_t inod) const; @@ -175,12 +181,6 @@ protected: // Internal utility methods // ======================== - //! \brief Returns a matrix with nodal coordinates for an element. - //! \param[in] iel Element index - //! \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; - //! \brief Returns the number of nodal points in each parameter direction. //! \param[out] n1 Number of nodes in first (u) direction //! \param[out] n2 Number of nodes in second (v) direction diff --git a/src/ASM/ASMs3D.C b/src/ASM/ASMs3D.C index 5863abd2..093cd3f5 100644 --- a/src/ASM/ASMs3D.C +++ b/src/ASM/ASMs3D.C @@ -75,6 +75,9 @@ bool ASMs3D::read (std::istream& is) return false; } + // Number of DOFs pr element + neldof = svol->order(0)*svol->order(1)*svol->order(2)*nf; + geo = svol; return true; } @@ -1097,6 +1100,25 @@ static double getElmSize (int p1, int p2, int p3, const Matrix& X) } +/*! + \brief Computes the characteristic element length from inverse Jacobian. +*/ + +static bool getG (const Matrix& Jinv, const Vector& du, Matrix& G) +{ + const size_t nsd = Jinv.cols(); + + G.resize(nsd,nsd,true); + + for (size_t k = 1;k <= nsd;k++) + for (size_t l = 1;l <= nsd;l++) + for (size_t m = 1;m <= nsd;m++) + G(k,l) += Jinv(m,k)*Jinv(m,l)/(du(k)*du(l)); + + return true; +} + + void ASMs3D::extractBasis (const Go::BasisDerivs& spline, Vector& N, Matrix& dNdu) { @@ -1335,10 +1357,19 @@ bool ASMs3D::integrate (Integrand& integrand, if (fe.detJxW == 0.0) continue; // skip singular points // Compute Hessian of coordinate mapping and 2nd order derivatives - if (integrand.getIntegrandType() == 2) + if (integrand.getIntegrandType() == 2) { if (!utl::Hessian(Hess,fe.d2NdX2,Jac,Xnod,d2Ndu2,dNdu)) return false; + // RUNAR + // Vector dXidu(3); +// dXidu(1) = this->getParametricLength(iel,1); +// dXidu(2) = this->getParametricLength(iel,2); +// dXidu(3) = this->getParametricLength(iel,3); +// if (!getG(Jac,dXidu,fe.G)) +// return false; + } + #if SP_DEBUG > 4 std::cout <<"\niel, ip = "<< iel <<" "<< ip <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX << std::endl; diff --git a/src/ASM/ASMs3D.h b/src/ASM/ASMs3D.h index 15f4de95..0b1806a0 100644 --- a/src/ASM/ASMs3D.h +++ b/src/ASM/ASMs3D.h @@ -118,6 +118,17 @@ public: //! This is used to reinitialize the patch after it has been refined. virtual void clear(bool retainGeometry = false); + //! \brief Returns a matrix with nodal coordinates for an element. + //! \param[in] iel Element index + //! \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; + + //! \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 + virtual void getNodalCoordinates(Matrix& X) const; + //! \brief Returns the global coordinates for the given node. //! \param[in] inod 1-based node index local to current patch virtual Vec3 getCoord(size_t inod) const; @@ -390,16 +401,6 @@ protected: //! \param[in] dir Local face index of the boundary face double getParametricArea(int iel, int dir) const; - //! \brief Returns a matrix with nodal coordinates for an element. - //! \param[in] iel Element index - //! \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; - //! \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 - virtual void getNodalCoordinates(Matrix& X) const; - //! \brief Returns the number of nodal points in each parameter direction. //! \param[out] n1 Number of nodes in first (u) direction //! \param[out] n2 Number of nodes in second (v) direction diff --git a/src/ASM/ASMs3DLag.h b/src/ASM/ASMs3DLag.h index 8c147fc8..bb51377e 100644 --- a/src/ASM/ASMs3DLag.h +++ b/src/ASM/ASMs3DLag.h @@ -138,18 +138,12 @@ public: virtual bool evalSolution(Matrix& sField, const Integrand& integrand, const RealArray* gpar, bool regular = true) const; -protected: - - // Internal utility methods - // ======================== - //! \brief Returns a matrix with nodal coordinates for an element. //! \param[in] iel Element index //! \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 diff --git a/src/ASM/ASMs3Dmx.C b/src/ASM/ASMs3Dmx.C index 89c71083..8bfaf037 100644 --- a/src/ASM/ASMs3Dmx.C +++ b/src/ASM/ASMs3Dmx.C @@ -153,6 +153,17 @@ bool ASMs3Dmx::generateFEMTopology () Go::BsplineBasis b1 = svol->basis(0).extendedBasis(svol->order(0)+1); Go::BsplineBasis b2 = svol->basis(1).extendedBasis(svol->order(1)+1); Go::BsplineBasis b3 = svol->basis(2).extendedBasis(svol->order(2)+1); + /* To lower order and regularity this can be used instead + std::vector::const_iterator first = ++svol->basis(0).begin(); + std::vector::const_iterator last = --svol->basis(0).end(); + Go::BsplineBasis b1 = Go::BsplineBasis(svol->order(0)-1,first,last); + first = ++svol->basis(1).begin(); + last = --svol->basis(1).end(); + Go::BsplineBasis b2 = Go::BsplineBasis(svol->order(1)-1,first,last); + first = ++svol->basis(2).begin(); + last = --svol->basis(2).end(); + Go::BsplineBasis b3 = Go::BsplineBasis(svol->order(2)-1,first,last); + */ // Note: Currently this is implemented for non-rational splines only. // TODO: Ask the splines people how to fix this properly, that is, how @@ -367,6 +378,11 @@ bool ASMs3Dmx::generateFEMTopology () } } + // Number of DOFs pr element for each basis + neldof1 = basis1->order(0)*basis1->order(1)*basis1->order(2)*nf1; + neldof2 = basis2->order(0)*basis2->order(1)*basis2->order(2)*nf2; + neldof = neldof1+neldof2; + #ifdef SP_DEBUG std::cout <<"NEL = "<< MLGE.size() <<" NNOD = "<< MLGN.size() << std::endl; #endif diff --git a/src/ASM/ASMs3Dmx.h b/src/ASM/ASMs3Dmx.h index f36d992a..167ddc4d 100644 --- a/src/ASM/ASMs3Dmx.h +++ b/src/ASM/ASMs3Dmx.h @@ -55,6 +55,12 @@ public: //! \param[in] retainGeometry If \e true, the spline geometry is not cleared. //! This is used to reinitialize the patch after it has been refined. virtual void clear(bool retainGeometry = false); + + //! \brief Returns a matrix with nodal coordinates for an element. + //! \param[in] iel Element index + //! \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; //! \brief Returns the global coordinates for the given node. //! \param[in] inod 1-based node index local to current patch @@ -174,12 +180,6 @@ protected: // Internal utility methods // ======================== - //! \brief Returns a matrix with nodal coordinates for an element. - //! \param[in] iel Element index - //! \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; - //! \brief Returns the number of nodal points in each parameter direction. //! \param[out] n1 Number of nodes in first (u) direction //! \param[out] n2 Number of nodes in second (v) direction