diff --git a/src/ASM/ASMmxBase.C b/src/ASM/ASMmxBase.C index 22b55a0f..9f9eb1f4 100644 --- a/src/ASM/ASMmxBase.C +++ b/src/ASM/ASMmxBase.C @@ -150,39 +150,12 @@ ASMmxBase::SurfaceVec ASMmxBase::establishBases (Go::SplineSurface* surf, else if (ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE) { result.resize(3); - - // basis1 should be one degree higher than basis2 and C^p-1 continuous - int ndim = surf->dimension(); - Go::BsplineBasis a1 = surf->basis(0); - Go::BsplineBasis a2 = surf->basis(1); - Go::BsplineBasis b1 = surf->basis(0).extendedBasis(surf->order_u()+1); - Go::BsplineBasis b2 = surf->basis(1).extendedBasis(surf->order_v()+1); - - // Compute parameter values of the Greville points - size_t i; - RealArray u0(a1.numCoefs()), v0(a2.numCoefs()); - for (i = 0; i < u0.size(); i++) - u0[i] = a1.grevilleParameter(i); - for (i = 0; i < v0.size(); i++) - v0[i] = a2.grevilleParameter(i); - RealArray ug(b1.numCoefs()), vg(b2.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); - - // Evaluate the spline surface at all points - // Project the coordinates onto the new basis (the 2nd XYZ is dummy here) - RealArray XYZ0(ndim*ug.size()*v0.size()), XYZ1(ndim*u0.size()*vg.size()); - surf->gridEvaluator(XYZ0,ug,v0); - surf->gridEvaluator(XYZ1,u0,vg); - result[2].reset(new Go::SplineSurface(*surf)); - result[0].reset(Go::SurfaceInterpolator::regularInterpolation(b1,a2, - ug,v0,XYZ0,ndim, - false,XYZ0)); - result[1].reset(Go::SurfaceInterpolator::regularInterpolation(a1,b2, - u0,vg,XYZ1,ndim, - false,XYZ1)); + result[0].reset(ASMmxBase::adjustBasis(*surf,{SplineUtils::AdjustOp::Original, + SplineUtils::AdjustOp::Lower})); + result[1].reset(ASMmxBase::adjustBasis(*surf,{SplineUtils::AdjustOp::Lower, + SplineUtils::AdjustOp::Original})); + result[2].reset(ASMmxBase::adjustBasis(*surf,{SplineUtils::AdjustOp::Lower, + SplineUtils::AdjustOp::Lower})); itgBasis = 3; } else if (type == SUBGRID) { // basis1 should be one degree higher than basis2 and C^p-1 continuous @@ -241,49 +214,18 @@ ASMmxBase::VolumeVec ASMmxBase::establishBases (Go::SplineVolume* svol, else if (ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE) { result.resize(4); - - // basis1 should be one degree higher than basis2 and C^p-1 continuous - int ndim = svol->dimension(); - Go::BsplineBasis a1 = svol->basis(0); - Go::BsplineBasis a2 = svol->basis(1); - Go::BsplineBasis a3 = svol->basis(2); - 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); - - // Compute parameter values of the Greville points - size_t i; - RealArray u0(a1.numCoefs()), v0(a2.numCoefs()), w0(a3.numCoefs()); - for (i = 0; i < u0.size(); i++) - u0[i] = a1.grevilleParameter(i); - for (i = 0; i < v0.size(); i++) - v0[i] = a2.grevilleParameter(i); - for (i = 0; i < w0.size(); i++) - w0[i] = a3.grevilleParameter(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 surface at all points - // Project the coordinates onto the new basis (the 2nd XYZ is dummy here) - RealArray XYZ0(ndim*ug.size()*v0.size()*w0.size()), XYZ1(ndim*u0.size()*vg.size()*w0.size()), XYZ2(ndim*u0.size()*v0.size()*wg.size()); - svol->gridEvaluator(ug,v0,w0,XYZ0); - svol->gridEvaluator(u0,vg,w0,XYZ1); - svol->gridEvaluator(u0,v0,wg,XYZ2); - result[0].reset(Go::VolumeInterpolator::regularInterpolation(b1,a2,a3, - ug,v0,w0,XYZ0,ndim, - false,XYZ0)); - result[1].reset(Go::VolumeInterpolator::regularInterpolation(a1,b2,a3, - u0,vg,w0,XYZ1,ndim, - false,XYZ1)); - result[2].reset(Go::VolumeInterpolator::regularInterpolation(a1,a2,b3, - u0,v0,wg,XYZ2,ndim, - false,XYZ2)); - result[3].reset(new Go::SplineVolume(*svol)); + result[0].reset(ASMmxBase::adjustBasis(*svol,{SplineUtils::AdjustOp::Original, + SplineUtils::AdjustOp::Lower, + SplineUtils::AdjustOp::Lower})); + result[1].reset(ASMmxBase::adjustBasis(*svol,{SplineUtils::AdjustOp::Lower, + SplineUtils::AdjustOp::Original, + SplineUtils::AdjustOp::Lower})); + result[2].reset(ASMmxBase::adjustBasis(*svol,{SplineUtils::AdjustOp::Lower, + SplineUtils::AdjustOp::Lower, + SplineUtils::AdjustOp::Original})); + result[3].reset(ASMmxBase::adjustBasis(*svol,{SplineUtils::AdjustOp::Lower, + SplineUtils::AdjustOp::Lower, + SplineUtils::AdjustOp::Lower})); itgBasis = 4; } else if (type == SUBGRID) { // basis1 should be one degree higher than basis2 and C^p-1 continuous diff --git a/src/ASM/ASMs2Dmx.C b/src/ASM/ASMs2Dmx.C index 4904babc..3eb06e36 100644 --- a/src/ASM/ASMs2Dmx.C +++ b/src/ASM/ASMs2Dmx.C @@ -206,17 +206,24 @@ bool ASMs2Dmx::generateFEMTopology () if (!surf) return false; if (m_basis.empty()) { + if (ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE && + (surf->order_u() < 3 || surf->order_v() < 3)) { + std::cerr << "*** RT basis cannot use a linear geometry." << std::endl; + return false; + } + m_basis = ASMmxBase::establishBases(surf, ASMmxBase::Type); // we need to project on something that is not one of our bases if (!projB) { if (ASMmxBase::Type == ASMmxBase::REDUCED_CONT_RAISE_BASIS1 || - ASMmxBase::Type == ASMmxBase::REDUCED_CONT_RAISE_BASIS2 || - ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE) + ASMmxBase::Type == ASMmxBase::REDUCED_CONT_RAISE_BASIS2) projB = ASMmxBase::adjustBasis(*surf,{SplineUtils::AdjustOp::Raise, SplineUtils::AdjustOp::Raise}); else if (ASMmxBase::Type == ASMmxBase::SUBGRID) projB = m_basis.front().get(); + else if (ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE) + projB = new Go::SplineSurface(*surf); else // FULL_CONT_RAISE_BASISx projB = m_basis[2-itgBasis].get(); } diff --git a/src/ASM/ASMs3Dmx.C b/src/ASM/ASMs3Dmx.C index 372443ef..864ed416 100644 --- a/src/ASM/ASMs3Dmx.C +++ b/src/ASM/ASMs3Dmx.C @@ -204,18 +204,24 @@ bool ASMs3Dmx::generateFEMTopology () if (!svol) return false; if (m_basis.empty()) { + if (ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE && + (svol->order(0) < 3 || svol->order(1) < 3 || svol->order(2) < 3)) { + std::cerr << "*** RT basis cannot use a linear geometry." << std::endl; + return false; + } m_basis = ASMmxBase::establishBases(svol, ASMmxBase::Type); // we need to project on something that is not one of our bases if (!projB) { if (ASMmxBase::Type == ASMmxBase::REDUCED_CONT_RAISE_BASIS1 || - ASMmxBase::Type == ASMmxBase::REDUCED_CONT_RAISE_BASIS2 || - ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE) + ASMmxBase::Type == ASMmxBase::REDUCED_CONT_RAISE_BASIS2) projB = ASMmxBase::adjustBasis(*svol,{SplineUtils::AdjustOp::Raise, SplineUtils::AdjustOp::Raise, SplineUtils::AdjustOp::Raise}); else if (ASMmxBase::Type == ASMmxBase::SUBGRID) projB = m_basis.front().get(); + else if (ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE) + projB = new Go::SplineVolume(*svol); else // FULL_CONT_RAISE_BASISx projB = m_basis[2-itgBasis].get(); } diff --git a/src/ASM/LR/ASMu2Dmx.C b/src/ASM/LR/ASMu2Dmx.C index ba180b08..a8ecac49 100644 --- a/src/ASM/LR/ASMu2Dmx.C +++ b/src/ASM/LR/ASMu2Dmx.C @@ -205,22 +205,30 @@ bool ASMu2Dmx::generateFEMTopology () } if (m_basis.empty()) { + if (ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE && + (tensorspline->order_u() < 3 || tensorspline->order_v() < 3)) { + std::cerr << "*** RT basis cannot use a linear geometry." << std::endl; + return false; + } + SurfaceVec svec = ASMmxBase::establishBases(tensorspline, ASMmxBase::Type); for (size_t b = 0; b < svec.size(); b++) m_basis.push_back(createLR(*svec[b])); - std::unique_ptr otherBasis( - ASMmxBase::adjustBasis(*tensorspline,{SplineUtils::AdjustOp::Raise, - SplineUtils::AdjustOp::Raise})); + std::unique_ptr otherBasis; + if (ASMmxBase::Type != ASMmxBase::DIV_COMPATIBLE) + otherBasis.reset(ASMmxBase::adjustBasis(*tensorspline,{SplineUtils::AdjustOp::Raise, + SplineUtils::AdjustOp::Raise})); // we need to project on something that is not one of our bases if (!projB) { if (ASMmxBase::Type == ASMmxBase::REDUCED_CONT_RAISE_BASIS1 || - ASMmxBase::Type == ASMmxBase::REDUCED_CONT_RAISE_BASIS2 || - ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE) + ASMmxBase::Type == ASMmxBase::REDUCED_CONT_RAISE_BASIS2) projB = createLR(*otherBasis); else if (ASMmxBase::Type == ASMmxBase::SUBGRID) projB = m_basis.front(); + else if (ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE) + projB = createLR(*tensorspline); else // FULL_CONT_RAISE_BASISx projB = m_basis[2-ASMmxBase::itgBasis]; } diff --git a/src/ASM/LR/ASMu3Dmx.C b/src/ASM/LR/ASMu3Dmx.C index d008db83..1a9ffa46 100644 --- a/src/ASM/LR/ASMu3Dmx.C +++ b/src/ASM/LR/ASMu3Dmx.C @@ -189,23 +189,32 @@ bool ASMu3Dmx::generateFEMTopology () return true; if (m_basis.empty()) { + if (ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE && + (tensorspline->order(0) < 3 || + tensorspline->order(1) < 3 || + tensorspline->order(2) < 3)) { + std::cerr << "*** RT basis cannot use a linear geometry." << std::endl; + return false; + } VolumeVec vvec = ASMmxBase::establishBases(tensorspline, ASMmxBase::Type); for (size_t b = 0; b < vvec.size(); b++) m_basis.push_back(std::make_shared(vvec[b].get())); - std::unique_ptr otherBasis( - ASMmxBase::adjustBasis(*tensorspline,{SplineUtils::AdjustOp::Raise, - SplineUtils::AdjustOp::Raise, - SplineUtils::AdjustOp::Raise})); + std::unique_ptr otherBasis; + if (ASMmxBase::Type != ASMmxBase::DIV_COMPATIBLE) + otherBasis.reset(ASMmxBase::adjustBasis(*tensorspline,{SplineUtils::AdjustOp::Raise, + SplineUtils::AdjustOp::Raise, + SplineUtils::AdjustOp::Raise})); // we need to project on something that is not one of our bases if (!projB) { if (ASMmxBase::Type == ASMmxBase::REDUCED_CONT_RAISE_BASIS1 || - ASMmxBase::Type == ASMmxBase::REDUCED_CONT_RAISE_BASIS2 || - ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE) + ASMmxBase::Type == ASMmxBase::REDUCED_CONT_RAISE_BASIS2) projB = std::make_shared(otherBasis.get()); else if (ASMmxBase::Type == ASMmxBase::SUBGRID) projB = m_basis.front(); + else if (ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE) + projB = std::make_shared(tensorspline); else // FULL_CONT_RAISE_BASISx projB = m_basis[2-ASMmxBase::itgBasis]; } diff --git a/src/ASM/LR/Test/TestASMu2Dmx.C b/src/ASM/LR/Test/TestASMu2Dmx.C index ad8caaca..c805173b 100644 --- a/src/ASM/LR/Test/TestASMu2Dmx.C +++ b/src/ASM/LR/Test/TestASMu2Dmx.C @@ -233,6 +233,7 @@ TEST(TestASMu2Dmx, WriteRT) ASMmxBase::Type = ASMmxBase::DIV_COMPATIBLE; ASMbase::resetNumbering(); ASMmxuSquare pch1({1,1,1}); + pch1.raiseOrder(1,1); EXPECT_TRUE(pch1.generateFEMTopology()); std::stringstream str; diff --git a/src/ASM/LR/Test/TestASMu3Dmx.C b/src/ASM/LR/Test/TestASMu3Dmx.C index b3c90f98..192b12f6 100644 --- a/src/ASM/LR/Test/TestASMu3Dmx.C +++ b/src/ASM/LR/Test/TestASMu3Dmx.C @@ -481,6 +481,7 @@ TEST(TestASMu3Dmx, WriteRT) ASMmxBase::Type = ASMmxBase::DIV_COMPATIBLE; ASMbase::resetNumbering(); ASMmxuCube pch1({1,1,1}); + pch1.raiseOrder(1,1,1,false); EXPECT_TRUE(pch1.generateFEMTopology()); std::stringstream str; diff --git a/src/ASM/LR/Test/TestLRSplineFields.C b/src/ASM/LR/Test/TestLRSplineFields.C index 59ee0a60..0d715881 100644 --- a/src/ASM/LR/Test/TestLRSplineFields.C +++ b/src/ASM/LR/Test/TestLRSplineFields.C @@ -57,6 +57,7 @@ TEST(TestLRSplineFields, Value2Dmx) { ASMmxBase::Type = ASMmxBase::DIV_COMPATIBLE; ASMmxuSquare patch({1,1,1}); + patch.raiseOrder(1,1); EXPECT_TRUE(patch.generateFEMTopology()); // {x+y+x*y, x-y+x*y} @@ -225,6 +226,7 @@ TEST(TestLRSplineFields, Value3Dmx) { ASMmxBase::Type = ASMmxBase::DIV_COMPATIBLE; ASMmxuCube patch({1,1,1,1}); + patch.raiseOrder(1,1,1,false); ASSERT_TRUE(patch.generateFEMTopology()); // {x+y+z+x*y*z, x+y-z+x*y*z, x-y+z+x*y*z} diff --git a/src/ASM/Test/TestASMs2Dmx.C b/src/ASM/Test/TestASMs2Dmx.C index 68006111..eabcb3db 100644 --- a/src/ASM/Test/TestASMs2Dmx.C +++ b/src/ASM/Test/TestASMs2Dmx.C @@ -188,6 +188,7 @@ TEST(TestASMs2Dmx, WriteRT) ASMmxBase::Type = ASMmxBase::DIV_COMPATIBLE; ASMbase::resetNumbering(); ASMmxSquare pch1({1,1,1}); + pch1.raiseOrder(1,1); EXPECT_TRUE(pch1.generateFEMTopology()); std::stringstream str; diff --git a/src/ASM/Test/TestASMs3Dmx.C b/src/ASM/Test/TestASMs3Dmx.C index 51624c57..10c7dfbd 100644 --- a/src/ASM/Test/TestASMs3Dmx.C +++ b/src/ASM/Test/TestASMs3Dmx.C @@ -415,6 +415,7 @@ TEST(TestASMs3Dmx, WriteRT) ASMmxBase::Type = ASMmxBase::DIV_COMPATIBLE; ASMbase::resetNumbering(); ASMmxCube pch1({1,1,1}); + pch1.raiseOrder(1,1,1); EXPECT_TRUE(pch1.generateFEMTopology()); std::stringstream str; diff --git a/src/ASM/Test/TestSplineFields.C b/src/ASM/Test/TestSplineFields.C index 4b93b885..81102979 100644 --- a/src/ASM/Test/TestSplineFields.C +++ b/src/ASM/Test/TestSplineFields.C @@ -52,6 +52,7 @@ TEST(TestSplineFields, Value2Dmx) { ASMmxBase::Type = ASMmxBase::DIV_COMPATIBLE; ASMmxSquare patch({1,1,1}); + patch.raiseOrder(1,1); EXPECT_TRUE(patch.generateFEMTopology()); // {x+y+x*y, x-y+x*y} @@ -205,6 +206,7 @@ TEST(TestSplineFields, Value3Dmx) { ASMmxBase::Type = ASMmxBase::DIV_COMPATIBLE; ASMmxCube patch({1,1,1,1}); + patch.raiseOrder(1,1,1); EXPECT_TRUE(patch.generateFEMTopology()); // {x+y+z+x*y*z, x+y-z+x*y*z, x-y+z+x*y*z}