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

addied: 3D COMPATIBLE

Browse Source
This commit is contained in:
timovanopstal
2016-01-15 16:10:56 +01:00
committed by Knut Morten Okstad
parent 4afe27a5d4
commit 736614ef37
3 changed files with 71 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

51
src/ASM/ASMmxBase.C
View File

@@ -295,6 +295,7 @@ ASMmxBase::VolumeVec ASMmxBase::establishBases(Go::SplineVolume* svol,
ug,vg,wg,XYZ,ndim,
false,XYZ));
}
result[1].reset(new Go::SplineVolume(*svol));
}
else if (type == REDUCED_CONT_RAISE_BASIS1 || type == REDUCED_CONT_RAISE_BASIS2)
{
@@ -302,8 +303,56 @@ ASMmxBase::VolumeVec ASMmxBase::establishBases(Go::SplineVolume* svol,
// 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));
}
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));
geoBasis = 4;
}
result[1].reset(new Go::SplineVolume(*svol));
if (type == FULL_CONT_RAISE_BASIS2 || type == REDUCED_CONT_RAISE_BASIS2)
std::swap(result[0], result[1]);

2
src/ASM/ASMs3D.C
View File

@@ -1297,7 +1297,7 @@ bool ASMs3D::updateDirichlet (const std::map<int,RealFunc*>& func,
for (size_t i = 0; i < dirich.size(); i++)
{
// Project the function onto the spline surface basis
Go::SplineSurface* dsurf = 0;
Go::SplineSurface* dsurf = nullptr;
if ((fit = func.find(dirich[i].code)) != func.end())
dsurf = SplineUtils::project(dirich[i].surf,*fit->second,time);
else if ((vfit = vfunc.find(dirich[i].code)) != vfunc.end())

24
src/ASM/ASMs3Dmx.C
View File

@@ -428,6 +428,8 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
PROFILE2("ASMs3Dmx::integrate(I)");
bool useElmVtx = integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS;
// Get Gaussian quadrature points and weights
const double* xg = GaussQuadrature::getCoord(nGauss);
const double* wg = GaussQuadrature::getWeight(nGauss);
@@ -493,6 +495,9 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
break;
}
if (useElmVtx)
this->getElementCorners(i1-1,i2-1,i3-1,fe.XC);
// Initialize element quantities
LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel,false);
if (!integrand.initElement(MNPC[iel-1],elem_sizes,nb,*A))
@@ -571,6 +576,8 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
PROFILE2("ASMs3Dmx::integrate(B)");
bool useElmVtx = integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS;
std::map<char,ThreadGroups>::const_iterator tit;
if ((tit = threadGroupsFace.find(lIndex)) == threadGroupsFace.end())
{
@@ -671,6 +678,9 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
break;
}
if (useElmVtx)
this->getElementCorners(i1-1,i2-1,i3-1,fe.XC);
// Initialize element quantities
LocalIntegral* A = integrand.getLocalIntegral(elem_sizes,fe.iel,true);
if (!integrand.initElementBou(MNPC[iel-1],elem_sizes,nb,*A))
@@ -843,7 +853,7 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const Vector& locSol,
Xtmp.multiply(splinex[b][i].basisValues,Ztmp);
Ytmp.insert(Ytmp.end(),Ztmp.begin(),Ztmp.end());
}
comp += nc[b];
comp += nc[b]*nb[b];
}
sField.fillColumn(1+i,Ytmp);
}
@@ -896,20 +906,26 @@ bool ASMs3Dmx::evalSolution (Matrix& sField, const IntegrandBase& integrand,
// Fetch indices of the non-zero basis functions at this point
std::vector<IntVec> ip(m_basis.size());
IntVec ipa;
size_t ofs = 0;
for (size_t b = 0; b < m_basis.size(); ++b) {
scatterInd(m_basis[b]->numCoefs(0),m_basis[b]->numCoefs(1),m_basis[b]->numCoefs(2),
m_basis[b]->order(0),m_basis[b]->order(1),m_basis[b]->order(2),
splinex[b][i].left_idx,ip[b]);
// Fetch associated control point coordinates
if (b == (size_t)geoBasis-1)
utl::gather(ip[geoBasis-1], 3, Xnod, Xtmp);
for (auto& it : ip[b])
it += ofs;
ipa.insert(ipa.end(), ip[b].begin(), ip[b].end());
ofs += nb[b];
}
fe.u = splinex[0][i].param[0];
fe.v = splinex[0][i].param[1];
fe.w = splinex[0][i].param[2];
// Fetch associated control point coordinates
utl::gather(ip[geoBasis-1], 3, Xnod, Xtmp);
// Fetch basis function derivatives at current integration point
for (size_t b = 0; b < m_basis.size(); ++b)
SplineUtils::extractBasis(splinex[b][i],fe.basis(b+1),dNxdu[b]);