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

changed: virtualize transferGaussPtVars / transferCtrlPtVars and implement for 3D

Browse Source
This commit is contained in:
Arne Morten Kvarving
2017-11-30 10:29:55 +01:00
parent c23ba91585
commit 982eb0ce3e
7 changed files with 330 additions and 32 deletions
Download Patch File Download Diff File Expand all files Collapse all files

31
src/ASM/ASMunstruct.h
View File

@@ -172,6 +172,37 @@ public:
virtual void remapErrors(RealArray& errors, const RealArray& orig,
bool elemErrors = false) const = 0;
//! \brief Transfers Gauss point variables from old basis to this patch.
//! \param[in] oldBasis The LR-spline basis to transfer from
//! \param[in] oldVars Gauss point variables associated with \a oldBasis
//! \param[out] newVars Gauss point variables associated with this patch.
//! \param[in] nGauss Number of Gauss points along a knot-span
virtual bool transferGaussPtVars(const LR::LRSpline* oldBasis,
const RealArray& oldVars, RealArray& newVars,
int nGauss) const = 0;
//! \brief Transfers Gauss point variables from old basis to this patch.
//! \param[in] oldBasis The LR-spline basis to transfer from
//! \param[in] oldVars Gauss point variables associated with \a oldBasis
//! \param[out] newVars Gauss point variables associated with this patch.
//! \param[in] nGauss Number of Gauss points along a knot-span
virtual bool transferGaussPtVarsN(const LR::LRSpline* oldBasis,
const RealArray& oldVars, RealArray& newVars,
int nGauss) const = 0;
//! \brief Transfers control point variables from old basis to this patch.
//! \param[in] oldBasis The LR-spline basis to transfer from
//! \param[out] newVars Gauss point variables associated with this patch.
//! \param[in] nGauss Number of Gauss points along a knot-span
virtual bool transferCntrlPtVars(const LR::LRSpline* oldBasis,
RealArray& newVars, int nGauss) const = 0;
//! \brief Transfers control point variables from old basis to this patch.
//! \param[in] oldBasis The LR-spline basis to transfer from
//! \param[in] oldVars Control point variables associated with \a oldBasis
//! \param[out] newVars Gauss point variables associated with this patch.
//! \param[in] nGauss Number of Gauss points along a knot-span
bool transferCntrlPtVars(LR::LRSpline* oldBasis,
const RealArray& oldVars, RealArray& newVars,
int nGauss, int nf = 1) const;
protected:
LR::LRSpline* geo; //!< Pointer to the actual spline geometry object

10
src/ASM/LR/ASMLRSpline.C
View File

@@ -384,3 +384,13 @@ void ASMunstruct::Sort(int u, int v, int orient,
return false;
});
}
bool ASMunstruct::transferCntrlPtVars (LR::LRSpline* oldBasis,
const RealArray& oldVars, RealArray& newVars,
int nGauss, int nf) const
{
oldBasis->rebuildDimension(nf);
oldBasis->setControlPoints(const_cast<RealArray&>(oldVars));
return this->transferCntrlPtVars(oldBasis,newVars,nGauss);
}

19
src/ASM/LR/ASMu2D.C
View File

@@ -1992,11 +1992,12 @@ size_t ASMu2D::getNoNodes (int) const
}
bool ASMu2D::transferGaussPtVars (const LR::LRSplineSurface* oldBasis,
bool ASMu2D::transferGaussPtVars (const LR::LRSpline* old_basis,
const RealArray& oldVars, RealArray& newVars,
int nGauss) const
{
const LR::LRSplineSurface* newBasis = this->getBasis();
const LR::LRSplineSurface* oldBasis = static_cast<const LR::LRSplineSurface*>(old_basis);
size_t nGp = nGauss*nGauss;
newVars.resize(newBasis->nElements()*nGp);
@@ -2044,11 +2045,12 @@ bool ASMu2D::transferGaussPtVars (const LR::LRSplineSurface* oldBasis,
}
bool ASMu2D::transferGaussPtVarsN (const LR::LRSplineSurface* oldBasis,
bool ASMu2D::transferGaussPtVarsN (const LR::LRSpline* old_basis,
const RealArray& oldVars, RealArray& newVars,
int nGauss) const
{
const LR::LRSplineSurface* newBasis = this->getBasis();
const LR::LRSplineSurface* oldBasis = static_cast<const LR::LRSplineSurface*>(old_basis);
size_t nGP = nGauss*nGauss;
newVars.clear();
@@ -2110,20 +2112,11 @@ bool ASMu2D::transferGaussPtVarsN (const LR::LRSplineSurface* oldBasis,
}
bool ASMu2D::transferCntrlPtVars (LR::LRSplineSurface* oldBasis,
const RealArray& oldVars, RealArray& newVars,
int nGauss) const
{
oldBasis->rebuildDimension(1);
oldBasis->setControlPoints(const_cast<RealArray&>(oldVars));
return this->transferCntrlPtVars(oldBasis,newVars,nGauss);
}
bool ASMu2D::transferCntrlPtVars (const LR::LRSplineSurface* oldBasis,
bool ASMu2D::transferCntrlPtVars (const LR::LRSpline* old_basis,
RealArray& newVars, int nGauss) const
{
const LR::LRSplineSurface* newBasis = this->getBasis();
const LR::LRSplineSurface* oldBasis = static_cast<const LR::LRSplineSurface*>(old_basis);
newVars.clear();
newVars.reserve(newBasis->nElements()*nGauss*nGauss*oldBasis->dimension());

32
src/ASM/LR/ASMu2D.h
View File

@@ -424,35 +424,29 @@ public:
double time) const;
//! \brief Transfers Gauss point variables from old basis to this patch.
//! \param[in] oldBasis The LR-spline basis to transfer from
//! \param[in] old_basis The LR-spline basis to transfer from
//! \param[in] oldVars Gauss point variables associated with \a oldBasis
//! \param[out] newVars Gauss point variables associated with this patch.
//! \param[in] nGauss Number of Gauss points along a knot-span
bool transferGaussPtVars(const LR::LRSplineSurface* oldBasis,
const RealArray& oldVars, RealArray& newVars,
int nGauss) const;
virtual bool transferGaussPtVars(const LR::LRSpline* old_basis,
const RealArray& oldVars, RealArray& newVars,
int nGauss) const;
//! \brief Transfers Gauss point variables from old basis to this patch.
//! \param[in] oldBasis The LR-spline basis to transfer from
//! \param[in] old_basis The LR-spline basis to transfer from
//! \param[in] oldVars Gauss point variables associated with \a oldBasis
//! \param[out] newVars Gauss point variables associated with this patch.
//! \param[in] nGauss Number of Gauss points along a knot-span
bool transferGaussPtVarsN(const LR::LRSplineSurface* oldBasis,
const RealArray& oldVars, RealArray& newVars,
int nGauss) const;
virtual bool transferGaussPtVarsN(const LR::LRSpline* old_basis,
const RealArray& oldVars, RealArray& newVars,
int nGauss) const;
using ASMunstruct::transferCntrlPtVars;
//! \brief Transfers control point variables from old basis to this patch.
//! \param[in] oldBasis The LR-spline basis to transfer from
//! \param[in] oldVars Control point variables associated with \a oldBasis
//! \param[in] old_basis The LR-spline basis to transfer from
//! \param[out] newVars Gauss point variables associated with this patch.
//! \param[in] nGauss Number of Gauss points along a knot-span
bool transferCntrlPtVars(LR::LRSplineSurface* oldBasis,
const RealArray& oldVars, RealArray& newVars,
int nGauss) const;
//! \brief Transfers control point variables from old basis to this patch.
//! \param[in] oldBasis The LR-spline basis to transfer from
//! \param[out] newVars Gauss point variables associated with this patch.
//! \param[in] nGauss Number of Gauss points along a knot-span
bool transferCntrlPtVars(const LR::LRSplineSurface* oldBasis,
RealArray& newVars, int nGauss) const;
virtual bool transferCntrlPtVars(const LR::LRSpline* old_basis,
RealArray& newVars, int nGauss) const;
protected:

175
src/ASM/LR/ASMu3D.C
View File

@@ -22,6 +22,7 @@
#include "TimeDomain.h"
#include "FiniteElement.h"
#include "GlobalIntegral.h"
#include "LagrangeInterpolator.h"
#include "LocalIntegral.h"
#include "IntegrandBase.h"
#include "CoordinateMapping.h"
@@ -2234,3 +2235,177 @@ bool ASMu3D::evaluate (const FunctionBase* func, RealArray& vec,
vec = ctrlPvals;
return ok;
}
bool ASMu3D::transferGaussPtVars (const LR::LRSpline* old_basis,
const RealArray& oldVars, RealArray& newVars,
int nGauss) const
{
const LR::LRSplineVolume* newBasis = this->getBasis();
const LR::LRSplineVolume* oldBasis = static_cast<const LR::LRSplineVolume*>(old_basis);
size_t nGp = nGauss*nGauss*nGauss;
newVars.resize(newBasis->nElements()*nGp);
const double* xi = GaussQuadrature::getCoord(nGauss);
LagrangeInterpolator interp(Vector(xi,nGauss));
for (int iEl = 0; iEl < newBasis->nElements(); iEl++)
{
const LR::Element* newEl = newBasis->getElement(iEl);
double u_center = 0.5*(newEl->umin() + newEl->umax());
double v_center = 0.5*(newEl->vmin() + newEl->vmax());
double w_center = 0.5*(newEl->wmin() + newEl->wmax());
int iOld = oldBasis->getElementContaining(u_center,v_center,w_center);
if (iOld < 0)
{
std::cerr <<" *** ASMu3D: Failed to locate element "<< iEl
<<" of the old mesh in the new mesh."<< std::endl;
return false;
}
const LR::Element* oldEl = oldBasis->getElement(iOld);
std::array<Matrix,3> I;
for (int i = 0; i < 3; i++)
{
RealArray UGP;
LR::getGaussPointParameters(newBasis, UGP, i, nGauss, iEl+1, xi);
double pmin = oldEl->getParmin(i);
double pmax = oldEl->getParmax(i);
for (size_t j = 0; j < UGP.size(); j++)
UGP[j] = -1.0 + 2.0*(UGP[j]-pmin)/(pmax-pmin);
// lagrangian interpolation
I[i] = interp.get(UGP);
}
Matrix data(nGauss,nGauss*nGauss);
data.fill(&oldVars[nGp*iOld],nGp);
Matrix I1 = I[0]*data;
for (int i = 0; i < nGauss; ++i)
for (int j = 0; j < nGauss; ++j)
for (int k = 0; k < nGauss; ++k)
data(j + 1, i + k*(nGauss) + 1) = I1(i+1, 1 + j + k*(nGauss));
I1 = I[1]*data;
for (int i = 0; i < nGauss; ++i)
for (int j = 0; j < nGauss; ++j)
for (int k = 0; k < nGauss; ++k)
data(k + 1, i + j*(nGauss) + 1) = I1(j+1, 1 + i + k*(nGauss));
I1 = I[2]*data;
for (int i = 0; i < nGauss; ++i)
for (int j = 0; j < nGauss; ++j)
for (int k = 0; k < nGauss; ++k)
newVars[iEl*nGp+i+(j + k*nGauss)*nGauss] = I1(k+1, 1+i + j*nGauss);
}
return true;
}
bool ASMu3D::transferGaussPtVarsN (const LR::LRSpline* old_basis,
const RealArray& oldVars, RealArray& newVars,
int nGauss) const
{
const LR::LRSplineVolume* newBasis = this->getBasis();
const LR::LRSplineVolume* oldBasis = static_cast<const LR::LRSplineVolume*>(old_basis);
size_t nGP = nGauss*nGauss*nGauss;
newVars.clear();
newVars.reserve(nGP*newBasis->nElements());
struct Param{ double u,v,w; };
std::vector<Param> oGP(nGP);
const double* xi = GaussQuadrature::getCoord(nGauss);
for (int iEl = 0; iEl < newBasis->nElements(); iEl++)
{
const LR::Element* newEl = newBasis->getElement(iEl);
double u_center = 0.5*(newEl->umin() + newEl->umax());
double v_center = 0.5*(newEl->vmin() + newEl->vmax());
double w_center = 0.5*(newEl->wmin() + newEl->wmax());
int iOld = oldBasis->getElementContaining(u_center,v_center,w_center);
if (iOld < 0)
{
std::cerr <<" *** ASMu3D: Failed to locate element "<< iEl
<<" of the old mesh in the new mesh."<< std::endl;
return false;
}
const LR::Element* oldEl = oldBasis->getElement(iOld);
// find parameters of old gauss points
double umin = oldEl->umin();
double vmin = oldEl->vmin();
double wmin = oldEl->wmin();
double du = 0.5*(oldEl->umax() - umin);
double dv = 0.5*(oldEl->vmax() - vmin);
double dw = 0.5*(oldEl->wmax() - wmin);
size_t l = 0;
for (int k = 0; k < nGauss; ++k)
for (int j = 0; j < nGauss; ++j)
for (int i = 0; i < nGauss; ++i, ++l) {
oGP[l].u = umin + du * (xi[i] + 1.0);
oGP[l].v = vmin + dv * (xi[j] + 1.0);
oGP[l].w = wmin + dw * (xi[k] + 1.0);
}
// parameters of new gauss points
umin = newEl->umin();
vmin = newEl->vmin();
wmin = newEl->wmin();
du = 0.5*(newEl->umax() - umin);
dv = 0.5*(newEl->vmax() - vmin);
dw = 0.5*(newEl->wmax() - wmin);
for (int k = 0; k < nGauss; ++k)
for (int j = 0; j < nGauss; ++j)
for (int i = 0; i < nGauss; ++i) {
double u = umin + du * (xi[i] + 1.0);
double v = vmin + dv * (xi[j] + 1.0);
double w = wmin + dw * (xi[k] + 1.0);
double dist = 1.0e16;
size_t near = 0;
for (size_t l = 0; l < nGP; ++l) {
Vec3 d(oGP[l].u-u, oGP[l].v-v, oGP[l].w-w);
double nd = d.length();
if (nd < dist) {
near = l;
dist = nd;
}
}
newVars.push_back(oldVars[iOld*nGP+near]);
}
}
return true;
}
bool ASMu3D::transferCntrlPtVars (const LR::LRSpline* old_basis,
RealArray& newVars, int nGauss) const
{
const LR::LRSplineVolume* newBasis = this->getBasis();
const LR::LRSplineVolume* oldBasis = static_cast<const LR::LRSplineVolume*>(old_basis);
newVars.clear();
newVars.reserve(newBasis->nElements()*nGauss*nGauss*nGauss*oldBasis->dimension());
const double* xi = GaussQuadrature::getCoord(nGauss);
for (int iEl = 0; iEl < newBasis->nElements(); iEl++)
{
RealArray U, V, W, ptVar;
LR::getGaussPointParameters(newBasis, U, 0, nGauss, iEl+1, xi);
LR::getGaussPointParameters(newBasis, V, 1, nGauss, iEl+1, xi);
LR::getGaussPointParameters(newBasis, W, 2, nGauss, iEl+1, xi);
for (int k = 0; k < nGauss; k++)
for (int j = 0; j < nGauss; j++)
for (int i = 0; i < nGauss; i++)
{
oldBasis->point(ptVar,U[i],V[j],W[k]);
for (size_t l = 0; l < ptVar.size(); l++)
newVars.push_back(ptVar[l]);
}
}
return true;
}

24
src/ASM/LR/ASMu3D.h
View File

@@ -339,6 +339,30 @@ public:
virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand,
const int* npe, char project = '\0') const;
//! \brief Transfers Gauss point variables from old basis to this patch.
//! \param[in] old_basis The LR-spline basis to transfer from
//! \param[in] oldVars Gauss point variables associated with \a oldBasis
//! \param[out] newVars Gauss point variables associated with this patch.
//! \param[in] nGauss Number of Gauss points along a knot-span
virtual bool transferGaussPtVars(const LR::LRSpline* old_basis,
const RealArray& oldVars, RealArray& newVars,
int nGauss) const;
//! \brief Transfers Gauss point variables from old basis to this patch.
//! \param[in] old_basis The LR-spline basis to transfer from
//! \param[in] oldVars Gauss point variables associated with \a oldBasis
//! \param[out] newVars Gauss point variables associated with this patch.
//! \param[in] nGauss Number of Gauss points along a knot-span
virtual bool transferGaussPtVarsN(const LR::LRSpline* old_basis,
const RealArray& oldVars, RealArray& newVars,
int nGauss) const;
using ASMunstruct::transferCntrlPtVars;
//! \brief Transfers control point variables from old basis to this patch.
//! \param[in] old_basis The LR-spline basis to transfer from
//! \param[out] newVars Gauss point variables associated with this patch.
//! \param[in] nGauss Number of Gauss points along a knot-span
virtual bool transferCntrlPtVars(const LR::LRSpline* old_basis,
RealArray& newVars, int nGauss) const;
private:
//! \brief Struct representing an inhomogeneous Dirichlet boundary condition.
struct DirichletFace

71
src/ASM/LR/Test/TestASMu3D.C
View File

@@ -11,7 +11,9 @@
//==============================================================================
#include "ASMu3D.h"
#include "GaussQuadrature.h"
#include "SIM3D.h"
#include "LRSpline/LRSplineVolume.h"
#include "gtest/gtest.h"
@@ -59,3 +61,72 @@ const std::vector<int> tests = {1,2,3,4,5,6};
INSTANTIATE_TEST_CASE_P(TestASMu3D,
TestASMu3D,
testing::ValuesIn(tests));
TEST(TestASMu3D, TransferGaussPtVars)
{
SIM3D sim(1);
sim.opt.discretization = ASM::LRSpline;
sim.createDefaultModel();
ASMu3D* pch = static_cast<ASMu3D*>(sim.getPatch(1));
size_t id[3];
const double* xi = GaussQuadrature::getCoord(3);
RealArray oldAr(3*3*3), newAr;
for (size_t idx = 0; idx < 3; ++idx) {
SIM3D simNew(1);
simNew.opt.discretization = ASM::LRSpline;
simNew.createDefaultModel();
ASMu3D* pchNew = static_cast<ASMu3D*>(simNew.getPatch(1));
pchNew->uniformRefine(idx, 1);
for (id[2] = 0; id[2] < 3; ++id[2])
for (id[1] = 0; id[1] < 3; ++id[1])
for (id[0] = 0; id[0] < 3; ++id[0])
oldAr[id[0]+(id[1]+id[2]*3)*3] = (1.0 + xi[id[idx]]) / 2.0;
pchNew->transferGaussPtVars(pch->getVolume(), oldAr, newAr, 3);
size_t k = 0;
for (size_t iEl = 0; iEl < 2; ++iEl)
for (id[2] = 0; id[2] < 3; ++id[2])
for (id[1] = 0; id[1] < 3; ++id[1])
for (id[0] = 0; id[0] < 3; ++id[0], ++k)
EXPECT_FLOAT_EQ(newAr[k], 0.5*iEl + 0.5*(xi[id[idx]] + 1.0) / 2.0);
}
}
TEST(TestASMu3D, TransferGaussPtVarsN)
{
SIM3D sim(1), sim2(1);
sim.opt.discretization = sim2.opt.discretization = ASM::LRSpline;
sim.createDefaultModel();
sim2.createDefaultModel();
ASMu3D* pch = static_cast<ASMu3D*>(sim.getPatch(1));
ASMu3D* pchNew = static_cast<ASMu3D*>(sim2.getPatch(1));
pchNew->uniformRefine(0, 1);
RealArray oldAr(3*3*3), newAr;
std::iota(oldAr.begin(), oldAr.end(), 1);
pchNew->transferGaussPtVarsN(pch->getVolume(), oldAr, newAr, 3);
static RealArray refAr = {{ 1.0, 1.0, 2.0,
4.0, 4.0, 5.0,
7.0, 7.0, 8.0,
10.0, 10.0, 11.0,
13.0, 13.0, 14.0,
16.0, 16.0, 17.0,
19.0, 19.0, 20.0,
22.0, 22.0, 23.0,
25.0, 25.0, 26.0,
2.0, 3.0, 3.0,
5.0, 6.0, 6.0,
8.0, 9.0, 9.0,
11.0, 12.0, 12.0,
14.0, 15.0, 15.0,
17.0, 18.0, 18.0,
20.0, 21.0, 21.0,
23.0, 24.0, 24.0,
26.0, 27.0, 27.0}};
EXPECT_EQ(refAr.size(), newAr.size());
for (size_t i = 0; i < refAr.size(); ++i)
EXPECT_FLOAT_EQ(refAr[i], newAr[i]);
}