diff --git a/src/ASM/ASMmxBase.C b/src/ASM/ASMmxBase.C index cbafbd70..f6e4bf33 100644 --- a/src/ASM/ASMmxBase.C +++ b/src/ASM/ASMmxBase.C @@ -14,6 +14,8 @@ #include "ASMmxBase.h" #include "GoTools/geometry/SplineSurface.h" #include "GoTools/geometry/SurfaceInterpolator.h" +#include "GoTools/trivariate/SplineVolume.h" +#include "GoTools/trivariate/VolumeInterpolator.h" char ASMmxBase::geoBasis = 2; ASMmxBase::MixedType ASMmxBase::Type = ASMmxBase::FULL_CONT_RAISE_BASIS1; @@ -184,3 +186,59 @@ ASMmxBase::SurfaceVec ASMmxBase::establishBases(Go::SplineSurface* surf, return result; } + + +ASMmxBase::VolumeVec ASMmxBase::establishBases(Go::SplineVolume* svol, + MixedType type) +{ + VolumeVec result(2); + // With mixed methods we need two separate spline spaces + if (type == FULL_CONT_RAISE_BASIS1 || type == FULL_CONT_RAISE_BASIS2) + { + // basis1 should be one degree higher than basis2 and C^p-1 continuous + int ndim = svol->dimension(); + 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 = ++surf->basis(0).begin(); + std::vector::const_iterator last = --surf->basis(0).end(); + Go::BsplineBasis b1 = Go::BsplineBasis(surf->order_u()-1,first,last); + first = ++surf->basis(1).begin(); + last = --surf->basis(1).end(); + Go::BsplineBasis b2 = Go::BsplineBasis(surf->order_v()-1,first,last); + */ + + // Compute parameter values of the Greville points + size_t i; + RealArray ug(b1.numCoefs()), vg(b2.numCoefs()), wg(b3.numCoefs()); + for (i = 0; i < ug.size(); i++) + ug[i] = b1.grevilleParameter(i); + for (i = 0; i < vg.size(); i++) + vg[i] = b2.grevilleParameter(i); + for (i = 0; i < wg.size(); i++) + wg[i] = b3.grevilleParameter(i); + + RealArray XYZ(ndim*ug.size()*vg.size()*wg.size()); + // Evaluate the spline surface at all points + svol->gridEvaluator(XYZ,ug,vg,wg); + + // Project the coordinates onto the new basis (the 2nd XYZ is dummy here) + result[0].reset(Go::VolumeInterpolator::regularInterpolation(b1,b2,b3, + ug,vg,wg,XYZ,ndim, + false,XYZ)); + } + else if (type == REDUCED_CONT_RAISE_BASIS1 || type == REDUCED_CONT_RAISE_BASIS2) + { + // Order-elevate basis1 such that it is of one degree higher than basis2 + // but only C^p-2 continuous + result[0].reset(new Go::SplineVolume(*svol)); + result[0]->raiseOrder(1,1,1); + } + result[1].reset(new Go::SplineVolume(*svol)); + + if (type == FULL_CONT_RAISE_BASIS2 || type == REDUCED_CONT_RAISE_BASIS2) + std::swap(result[0], result[1]); + + return result; +} diff --git a/src/ASM/ASMmxBase.h b/src/ASM/ASMmxBase.h index 83d75d1b..8f736e5a 100644 --- a/src/ASM/ASMmxBase.h +++ b/src/ASM/ASMmxBase.h @@ -19,6 +19,7 @@ namespace Go { class SplineSurface; + class SplineVolume; } @@ -74,6 +75,7 @@ public: static MixedType Type; //!< Type of mixed formulation used typedef std::vector> SurfaceVec; //!< Convenience type + typedef std::vector> VolumeVec; //!< Convenience type //! \brief Establish mixed bases //! \param[in] surf The base basis to use. @@ -81,6 +83,12 @@ public: //! \return Vector with bases. static SurfaceVec establishBases(Go::SplineSurface* surf, MixedType type); + //! \brief Establish mixed bases + //! \param[in] vol The base basis to use. + //! \param[in] type The type of bases to establish. + //! \return Vector with bases. + static VolumeVec establishBases(Go::SplineVolume* vol, MixedType type); + private: std::vector MADOF; //!< Matrix of accumulated DOFs for this patch diff --git a/src/ASM/ASMs3Dmx.C b/src/ASM/ASMs3Dmx.C index e3e66b40..88508403 100644 --- a/src/ASM/ASMs3Dmx.C +++ b/src/ASM/ASMs3Dmx.C @@ -175,73 +175,7 @@ bool ASMs3Dmx::generateFEMTopology () if (!svol) return false; if (m_basis.empty()) - { - m_basis.resize(2); - // With mixed methods we need two separate spline spaces - if (ASMmxBase::Type == FULL_CONT_RAISE_BASIS1 || - ASMmxBase::Type == FULL_CONT_RAISE_BASIS2) - { - // basis1 should be one degree higher than basis2 and C^p-1 continuous - int ndim = svol->dimension(); - 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 - // may we obtain the correct weights for basis1 when *svol is a NURBS? - if (svol->rational()) - std::cout <<"WARNING: The geometry basis is rational (using NURBS).\n" - <<" The basis for the unknown fields of one degree" - <<" higher will however be non-rational.\n" - <<" This may affect accuracy.\n"<< std::endl; - - // Compute parameter values of the Greville points of the new basis - size_t i; - RealArray ug(b1.numCoefs()), vg(b2.numCoefs()), wg(b3.numCoefs()); - for (i = 0; i < ug.size(); i++) - ug[i] = b1.grevilleParameter(i); - for (i = 0; i < vg.size(); i++) - vg[i] = b2.grevilleParameter(i); - for (i = 0; i < wg.size(); i++) - wg[i] = b3.grevilleParameter(i); - - // Evaluate the spline volume at all points - RealArray XYZ(ndim*ug.size()*vg.size()*wg.size()); - svol->gridEvaluator(ug,vg,wg,XYZ); - - // Project the coordinates onto the new basis (the 2nd XYZ is dummy here) - m_basis[0].reset(Go::VolumeInterpolator::regularInterpolation(b1,b2,b3, - ug,vg,wg,XYZ,ndim, - false,XYZ)); - } - else if (ASMmxBase::Type == REDUCED_CONT_RAISE_BASIS1 || - ASMmxBase::Type == REDUCED_CONT_RAISE_BASIS2) - { - // Order-elevate basis1 such that it is of one degree higher than basis2 - // but only C^p-2 continuous - m_basis[0].reset(new Go::SplineVolume(*svol)); - m_basis[0]->raiseOrder(1,1,1); - } - m_basis[1].reset(new Go::SplineVolume(*svol)); - - if (ASMmxBase::Type == FULL_CONT_RAISE_BASIS2 || - ASMmxBase::Type == REDUCED_CONT_RAISE_BASIS2) { - std::swap(m_basis[0], m_basis[1]); - svol = m_basis[0].get(); - } - } + m_basis = ASMmxBase::establishBases(svol, ASMmxBase::Type); nb.clear(); for (auto it : m_basis)