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

changed: establishing RT basis now uses order lowering

Browse Source 
does not make sense to use order raising since the geometry
needs to be the highest order
This commit is contained in:
Arne Morten Kvarving
2023-09-07 23:00:04 +02:00
parent eaefc70d26
commit dca21c8856
11 changed files with 71 additions and 91 deletions
Download Patch File Download Diff File Expand all files Collapse all files

94
src/ASM/ASMmxBase.C
View File

@@ -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

11
src/ASM/ASMs2Dmx.C
View File

@@ -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();
}

10
src/ASM/ASMs3Dmx.C
View File

@@ -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();
}

18
src/ASM/LR/ASMu2Dmx.C
View File

@@ -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<Go::SplineSurface> otherBasis(
ASMmxBase::adjustBasis(*tensorspline,{SplineUtils::AdjustOp::Raise,
SplineUtils::AdjustOp::Raise}));
std::unique_ptr<Go::SplineSurface> 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];
}

21
src/ASM/LR/ASMu3Dmx.C
View File

@@ -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<LR::LRSplineVolume>(vvec[b].get()));
std::unique_ptr<Go::SplineVolume> otherBasis(
ASMmxBase::adjustBasis(*tensorspline,{SplineUtils::AdjustOp::Raise,
SplineUtils::AdjustOp::Raise,
SplineUtils::AdjustOp::Raise}));
std::unique_ptr<Go::SplineVolume> 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<LR::LRSplineVolume>(otherBasis.get());
else if (ASMmxBase::Type == ASMmxBase::SUBGRID)
projB = m_basis.front();
else if (ASMmxBase::Type == ASMmxBase::DIV_COMPATIBLE)
projB = std::make_shared<LR::LRSplineVolume>(tensorspline);
else // FULL_CONT_RAISE_BASISx
projB = m_basis[2-ASMmxBase::itgBasis];
}

1
src/ASM/LR/Test/TestASMu2Dmx.C
View File

@@ -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;

1
src/ASM/LR/Test/TestASMu3Dmx.C
View File

@@ -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;

2
src/ASM/LR/Test/TestLRSplineFields.C
View File

@@ -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}

1
src/ASM/Test/TestASMs2Dmx.C
View File

@@ -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;

1
src/ASM/Test/TestASMs3Dmx.C
View File

@@ -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;

2
src/ASM/Test/TestSplineFields.C
View File

@@ -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}