added: support for using fields as anasol functions
This commit is contained in:
Arne Morten Kvarving
parent
7c7121b25c
commit
59cefa93b1
@ -61,9 +61,16 @@ macro(IFEM_add_unittests IFEM_PATH)
|
||||
list(REMOVE_ITEM TEST_SOURCES ${IFEM_PATH}/src/LinAlg/Test/TestISTLPETScMatrix.C)
|
||||
endif()
|
||||
|
||||
if(NOT HDF5_FOUND)
|
||||
list(REMOVE_ITEM TEST_SOURCES ${IFEM_PATH}/src/Utility/Test/TestFieldFunctions.C)
|
||||
list(REMOVE_ITEM TEST_SOURCES ${IFEM_PATH}/src/Utility/Test/TestFieldFunctionsLR.C)
|
||||
endif()
|
||||
|
||||
if(LRSPLINE_FOUND OR LRSpline_FOUND)
|
||||
file(GLOB LR_TEST_SRCS ${PROJECT_SOURCE_DIR}/src/ASM/LR/Test/*.C)
|
||||
list(APPEND TEST_SOURCES ${LR_TEST_SRCS})
|
||||
else()
|
||||
list(REMOVE_ITEM TEST_SOURCES ${IFEM_PATH}/src/Utility/Test/TestFieldFunctionsLR.C)
|
||||
endif()
|
||||
|
||||
IFEM_add_test_app("${TEST_SOURCES}"
|
||||
|
||||
@ -783,7 +783,8 @@ bool ASMs1D::integrate (Integrand& integrand,
|
||||
FiniteElement fe(p1);
|
||||
Matrix dNdu, Jac;
|
||||
Matrix3D d2Ndu2, Hess;
|
||||
Vec4 X;
|
||||
double param[3] = { 0.0, 0.0, 0.0 };
|
||||
Vec4 X(param);
|
||||
|
||||
if (nsd > 1 && (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES))
|
||||
fe.G.resize(nsd,2); // For storing d{X}/du and d2{X}/du2
|
||||
@ -825,8 +826,8 @@ bool ASMs1D::integrate (Integrand& integrand,
|
||||
// Local element coordinates of current integration point
|
||||
fe.xi = xr[i];
|
||||
|
||||
// Parameter values of current integration point
|
||||
fe.u = redpar(1+i,1+iel);
|
||||
// Parameter values of current integration point
|
||||
fe.u = param[0] = redpar(1+i,1+iel);
|
||||
|
||||
if (integrand.getIntegrandType() & Integrand::NO_DERIVATIVES)
|
||||
this->extractBasis(fe.u,fe.N);
|
||||
@ -840,7 +841,7 @@ bool ASMs1D::integrate (Integrand& integrand,
|
||||
}
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = fe.Xn * fe.N;
|
||||
X.assign(fe.Xn * fe.N);
|
||||
X.t = time.t;
|
||||
|
||||
// Compute the reduced integration terms of the integrand
|
||||
@ -860,7 +861,7 @@ bool ASMs1D::integrate (Integrand& integrand,
|
||||
fe.xi = xg[i];
|
||||
|
||||
// Parameter value of current integration point
|
||||
fe.u = gpar(1+i,1+iel);
|
||||
fe.u = param[0] = gpar(1+i,1+iel);
|
||||
|
||||
// Compute basis functions and derivatives
|
||||
if (integrand.getIntegrandType() & Integrand::NO_DERIVATIVES)
|
||||
@ -894,7 +895,7 @@ bool ASMs1D::integrate (Integrand& integrand,
|
||||
}
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = fe.Xn * fe.N;
|
||||
X.assign(fe.Xn * fe.N);
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
|
||||
@ -1547,7 +1547,8 @@ bool ASMs2D::integrate (Integrand& integrand,
|
||||
Matrix dNdu, Xnod, Jac;
|
||||
Matrix3D d2Ndu2, Hess;
|
||||
double dXidu[2];
|
||||
Vec4 X;
|
||||
double param[3] = { 0.0, 0.0, 0.0 };
|
||||
Vec4 X(param);
|
||||
for (size_t i = 0; i < threadGroups[g][t].size() && ok; i++)
|
||||
{
|
||||
int iel = threadGroups[g][t][i];
|
||||
@ -1652,8 +1653,8 @@ bool ASMs2D::integrate (Integrand& integrand,
|
||||
fe.eta = xr[j];
|
||||
|
||||
// Parameter values of current integration point
|
||||
fe.u = redpar[0](i+1,i1-p1+1);
|
||||
fe.v = redpar[1](j+1,i2-p2+1);
|
||||
fe.u = param[0] = redpar[0](i+1,i1-p1+1);
|
||||
fe.v = param[1] = redpar[1](j+1,i2-p2+1);
|
||||
|
||||
// Fetch basis function derivatives at current point
|
||||
SplineUtils::extractBasis(splineRed[ip],fe.N,dNdu);
|
||||
@ -1663,7 +1664,7 @@ bool ASMs2D::integrate (Integrand& integrand,
|
||||
if (fe.detJxW == 0.0) continue; // skip singular points
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.N;
|
||||
X.assign(Xnod * fe.N);
|
||||
X.t = time.t;
|
||||
|
||||
// Compute the reduced integration terms of the integrand
|
||||
@ -1688,8 +1689,8 @@ bool ASMs2D::integrate (Integrand& integrand,
|
||||
fe.eta = xg[j];
|
||||
|
||||
// Parameter values of current integration point
|
||||
fe.u = gpar[0](i+1,i1-p1+1);
|
||||
fe.v = gpar[1](j+1,i2-p2+1);
|
||||
fe.u = param[0] = gpar[0](i+1,i1-p1+1);
|
||||
fe.v = param[1] = gpar[1](j+1,i2-p2+1);
|
||||
|
||||
// Fetch basis function derivatives at current integration point
|
||||
if (use2ndDer)
|
||||
@ -1721,7 +1722,7 @@ bool ASMs2D::integrate (Integrand& integrand,
|
||||
#endif
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.N;
|
||||
X.assign(Xnod * fe.N);
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
@ -1847,7 +1848,7 @@ bool ASMs2D::integrate (Integrand& integrand,
|
||||
{
|
||||
// Compute the element center
|
||||
fe.h = this->getElementCorners(i1-1,i2-1,fe.XC);
|
||||
X = 0.25*(fe.XC[0]+fe.XC[1]+fe.XC[2]+fe.XC[3]);
|
||||
X.assign(0.25*(fe.XC[0]+fe.XC[1]+fe.XC[2]+fe.XC[3]));
|
||||
}
|
||||
else if (useElmVtx)
|
||||
fe.h = this->getElementCorners(i1-1,i2-1,fe.XC);
|
||||
@ -1908,6 +1909,7 @@ bool ASMs2D::integrate (Integrand& integrand,
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.N;
|
||||
X.u = elmPts[ip].data();
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
@ -1963,7 +1965,8 @@ bool ASMs2D::integrate (Integrand& integrand,
|
||||
FiniteElement fe(p1*p2);
|
||||
Matrix dNdu, Xnod, Jac;
|
||||
Vector dN;
|
||||
Vec4 X;
|
||||
double param[3] = { 0.0, 0.0, 0.0 };
|
||||
Vec4 X(param);
|
||||
Vec3 normal;
|
||||
double u[2], v[2];
|
||||
bool hasInterfaceElms = MLGE.size() > nel && MLGE.size() != 2*nel;
|
||||
@ -2031,16 +2034,16 @@ bool ASMs2D::integrate (Integrand& integrand,
|
||||
{
|
||||
fe.xi = edgeDir;
|
||||
fe.eta = xg[i];
|
||||
fe.u = edgeDir > 0 ? u[1] : u[0];
|
||||
fe.v = 0.5*((v[1]-v[0])*xg[i] + v[1]+v[0]);
|
||||
fe.u = param[0] = edgeDir > 0 ? u[1] : u[0];
|
||||
fe.v = param[1] = 0.5*((v[1]-v[0])*xg[i] + v[1]+v[0]);
|
||||
fe.p = p1 - 1;
|
||||
}
|
||||
else
|
||||
{
|
||||
fe.xi = xg[i];
|
||||
fe.eta = edgeDir/2;
|
||||
fe.u = 0.5*((u[1]-u[0])*xg[i] + u[1]+u[0]);
|
||||
fe.v = edgeDir > 0 ? v[1] : v[0];
|
||||
fe.u = param[0] = 0.5*((u[1]-u[0])*xg[i] + u[1]+u[0]);
|
||||
fe.v = param[1] = edgeDir > 0 ? v[1] : v[0];
|
||||
fe.p = p2 - 1;
|
||||
}
|
||||
|
||||
@ -2054,7 +2057,7 @@ bool ASMs2D::integrate (Integrand& integrand,
|
||||
if (edgeDir < 0) normal *= -1.0;
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.N;
|
||||
X.assign(Xnod * fe.N);
|
||||
X.t = time.t;
|
||||
|
||||
if (integrand.getIntegrandType() & Integrand::NORMAL_DERIVS)
|
||||
@ -2163,9 +2166,10 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex,
|
||||
fe.xi = fe.eta = edgeDir < 0 ? -1.0 : 1.0;
|
||||
fe.u = gpar[0](1,1);
|
||||
fe.v = gpar[1](1,1);
|
||||
double param[3] = { fe.u, fe.v, 0.0 };
|
||||
|
||||
Matrix dNdu, Xnod, Jac;
|
||||
Vec4 X;
|
||||
Vec4 X(param);
|
||||
Vec3 normal;
|
||||
double dXidu[2];
|
||||
|
||||
@ -2229,12 +2233,12 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex,
|
||||
if (gpar[0].size() > 1)
|
||||
{
|
||||
fe.xi = xg[i];
|
||||
fe.u = gpar[0](i+1,i1-p1+1);
|
||||
fe.u = param[0] = gpar[0](i+1,i1-p1+1);
|
||||
}
|
||||
if (gpar[1].size() > 1)
|
||||
{
|
||||
fe.eta = xg[i];
|
||||
fe.v = gpar[1](i+1,i2-p2+1);
|
||||
fe.v = param[1] = gpar[1](i+1,i2-p2+1);
|
||||
}
|
||||
|
||||
// Fetch basis function derivatives at current integration point
|
||||
@ -2251,7 +2255,7 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex,
|
||||
utl::getGmat(Jac,dXidu,fe.G);
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.N;
|
||||
X.assign(Xnod * fe.N);
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
@ -2348,9 +2352,11 @@ bool ASMs2D::tesselate (ElementBlock& grid, const int* npe) const
|
||||
// Establish the block grid coordinates
|
||||
size_t i, j, l;
|
||||
grid.resize(nx,ny);
|
||||
for (i = l = 0; i < grid.getNoNodes(); i++, l += surf->dimension())
|
||||
for (i = l = 0; i < grid.getNoNodes(); i++, l += surf->dimension()) {
|
||||
grid.setParams(i, gpar[0][i % nx], gpar[1][i / nx]);
|
||||
for (j = 0; j < nsd; j++)
|
||||
grid.setCoor(i,j,XYZ[l+j]);
|
||||
}
|
||||
|
||||
// Establish the block grid topology
|
||||
int ie, nse1 = npe[0] - 1;
|
||||
|
||||
@ -543,7 +543,8 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
|
||||
Matrix3D Hess;
|
||||
double dXidu[2];
|
||||
Matrix Xnod, Jac;
|
||||
Vec4 X;
|
||||
double param[3] = { 0.0, 0.0, 0.0 };
|
||||
Vec4 X(param);
|
||||
for (size_t i = 0; i < threadGroups[g][t].size() && ok; ++i)
|
||||
{
|
||||
int iel = threadGroups[g][t][i];
|
||||
@ -602,8 +603,8 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
|
||||
fe.eta = xg[j];
|
||||
|
||||
// Parameter values of current integration point
|
||||
fe.u = gpar[0](i+1,i1-p1+1);
|
||||
fe.v = gpar[1](j+1,i2-p2+1);
|
||||
fe.u = param[0] = gpar[0](i+1,i1-p1+1);
|
||||
fe.v = param[1] = gpar[1](j+1,i2-p2+1);
|
||||
|
||||
// Fetch basis function derivatives at current integration point
|
||||
for (size_t b = 0; b < m_basis.size(); ++b)
|
||||
@ -638,7 +639,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
|
||||
utl::getGmat(Jac,dXidu,fe.G);
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.basis(geoBasis);
|
||||
X.assign(Xnod * fe.basis(geoBasis));
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
@ -802,7 +803,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
|
||||
if (edgeDir < 0) normal *= -1.0;
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.basis(geoBasis);
|
||||
X .assign(Xnod * fe.basis(geoBasis));
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
|
||||
@ -1737,7 +1737,8 @@ bool ASMs3D::integrate (Integrand& integrand,
|
||||
Matrix dNdu, Xnod, Jac;
|
||||
Matrix3D d2Ndu2, Hess;
|
||||
double dXidu[3];
|
||||
Vec4 X;
|
||||
double param[3];
|
||||
Vec4 X(param);
|
||||
for (size_t l = 0; l < threadGroupsVol[g][t].size() && ok; l++)
|
||||
{
|
||||
int iel = threadGroupsVol[g][t][l];
|
||||
@ -1848,9 +1849,9 @@ bool ASMs3D::integrate (Integrand& integrand,
|
||||
fe.zeta = xr[k];
|
||||
|
||||
// Parameter values of current integration point
|
||||
fe.u = redpar[0](i+1,i1-p1+1);
|
||||
fe.v = redpar[1](j+1,i2-p2+1);
|
||||
fe.w = redpar[2](k+1,i3-p3+1);
|
||||
fe.u = param[0] = redpar[0](i+1,i1-p1+1);
|
||||
fe.v = param[1] = redpar[1](j+1,i2-p2+1);
|
||||
fe.w = param[2] = redpar[2](k+1,i3-p3+1);
|
||||
|
||||
// Fetch basis function derivatives at current point
|
||||
SplineUtils::extractBasis(splineRed[ip],fe.N,dNdu);
|
||||
@ -1886,9 +1887,9 @@ bool ASMs3D::integrate (Integrand& integrand,
|
||||
fe.zeta = xg[k];
|
||||
|
||||
// Parameter values of current integration point
|
||||
fe.u = gpar[0](i+1,i1-p1+1);
|
||||
fe.v = gpar[1](j+1,i2-p2+1);
|
||||
fe.w = gpar[2](k+1,i3-p3+1);
|
||||
fe.u = param[0] = gpar[0](i+1,i1-p1+1);
|
||||
fe.v = param[1] = gpar[1](j+1,i2-p2+1);
|
||||
fe.w = param[1] = gpar[2](k+1,i3-p3+1);
|
||||
|
||||
// Fetch basis function derivatives at current integration point
|
||||
if (use2ndDer)
|
||||
@ -1920,7 +1921,7 @@ bool ASMs3D::integrate (Integrand& integrand,
|
||||
#endif
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.N;
|
||||
X.assign(Xnod * fe.N);
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
@ -2114,6 +2115,7 @@ bool ASMs3D::integrate (Integrand& integrand,
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.N;
|
||||
X.u = itgPts[iel][ip].data();
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
@ -2254,7 +2256,8 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex,
|
||||
fe.w = gpar[2](1,1);
|
||||
|
||||
Matrix dNdu, Xnod, Jac;
|
||||
Vec4 X;
|
||||
double param[3] = { fe.u, fe.v, fe.w };
|
||||
Vec4 X(param);
|
||||
Vec3 normal;
|
||||
double dXidu[3];
|
||||
|
||||
@ -2341,17 +2344,17 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex,
|
||||
if (gpar[0].size() > 1)
|
||||
{
|
||||
fe.xi = xg[k1];
|
||||
fe.u = gpar[0](k1+1,i1-p1+1);
|
||||
fe.u = param[0] = gpar[0](k1+1,i1-p1+1);
|
||||
}
|
||||
if (gpar[1].size() > 1)
|
||||
{
|
||||
fe.eta = xg[k2];
|
||||
fe.v = gpar[1](k2+1,i2-p2+1);
|
||||
fe.v = param[1] = gpar[1](k2+1,i2-p2+1);
|
||||
}
|
||||
if (gpar[2].size() > 1)
|
||||
{
|
||||
fe.zeta = xg[k3];
|
||||
fe.w = gpar[2](k3+1,i3-p3+1);
|
||||
fe.w = param[2] = gpar[2](k3+1,i3-p3+1);
|
||||
}
|
||||
|
||||
// Fetch basis function derivatives at current integration point
|
||||
@ -2368,7 +2371,7 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex,
|
||||
utl::getGmat(Jac,dXidu,fe.G);
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.N;
|
||||
X.assign(Xnod * fe.N);
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
@ -2472,7 +2475,8 @@ bool ASMs3D::integrateEdge (Integrand& integrand, int lEdge,
|
||||
if (gpar[2].size() == 1) fe.zeta = fe.w == svol->startparam(2) ? -1.0 : 1.0;
|
||||
|
||||
Matrix dNdu, Xnod, Jac;
|
||||
Vec4 X;
|
||||
double param[3] = { 0.0, 0.0, 0.0 };
|
||||
Vec4 X(param);
|
||||
Vec3 tang;
|
||||
|
||||
|
||||
@ -2542,9 +2546,9 @@ bool ASMs3D::integrateEdge (Integrand& integrand, int lEdge,
|
||||
for (int i = 0; i < nGauss && ok; i++, ip++, fe.iGP++)
|
||||
{
|
||||
// Parameter values of current integration point
|
||||
if (gpar[0].size() > 1) fe.u = gpar[0](i+1,i1-p1+1);
|
||||
if (gpar[1].size() > 1) fe.v = gpar[1](i+1,i2-p2+1);
|
||||
if (gpar[2].size() > 1) fe.w = gpar[2](i+1,i3-p3+1);
|
||||
if (gpar[0].size() > 1) fe.u = param[0] = gpar[0](i+1,i1-p1+1);
|
||||
if (gpar[1].size() > 1) fe.v = param[1] = gpar[1](i+1,i2-p2+1);
|
||||
if (gpar[2].size() > 1) fe.w = param[2] = gpar[2](i+1,i3-p3+1);
|
||||
|
||||
// Fetch basis function derivatives at current integration point
|
||||
SplineUtils::extractBasis(spline[ip],fe.N,dNdu);
|
||||
@ -2554,7 +2558,7 @@ bool ASMs3D::integrateEdge (Integrand& integrand, int lEdge,
|
||||
if (fe.detJxW == 0.0) continue; // skip singular points
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.N;
|
||||
X.assign(Xnod * fe.N);
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
@ -2648,8 +2652,10 @@ bool ASMs3D::tesselate (ElementBlock& grid, const int* npe) const
|
||||
// Establish the block grid coordinates
|
||||
size_t i, j, k, l;
|
||||
grid.resize(nx,ny,nz);
|
||||
for (i = j = 0; i < grid.getNoNodes(); i++, j += svol->dimension())
|
||||
for (i = j = 0; i < grid.getNoNodes(); i++, j += svol->dimension()) {
|
||||
grid.setParams(i,gpar[0][i % nx], gpar[1][i/nx % ny], gpar[2][i / (nx*ny)]);
|
||||
grid.setCoor(i,XYZ[j],XYZ[j+1],XYZ[j+2]);
|
||||
}
|
||||
|
||||
// Establish the block grid topology
|
||||
int ie, nse1 = npe[0] - 1;
|
||||
|
||||
@ -489,7 +489,8 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
|
||||
Matrix3D Hess;
|
||||
double dXidu[3];
|
||||
Matrix Xnod, Jac;
|
||||
Vec4 X;
|
||||
double param[3];
|
||||
Vec4 X(param);
|
||||
for (size_t l = 0; l < threadGroupsVol[g][t].size() && ok; ++l)
|
||||
{
|
||||
int iel = threadGroupsVol[g][t][l];
|
||||
@ -552,9 +553,9 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
|
||||
fe.zeta = xg[k];
|
||||
|
||||
// Parameter values of current integration point
|
||||
fe.u = gpar[0](i+1,i1-p1+1);
|
||||
fe.v = gpar[1](j+1,i2-p2+1);
|
||||
fe.w = gpar[2](k+1,i3-p3+1);
|
||||
fe.u = param[0] = gpar[0](i+1,i1-p1+1);
|
||||
fe.v = param[1] = gpar[1](j+1,i2-p2+1);
|
||||
fe.w = param[2] = gpar[2](k+1,i3-p3+1);
|
||||
|
||||
// Fetch basis function derivatives at current integration point
|
||||
if (use2ndDer)
|
||||
@ -590,7 +591,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
|
||||
utl::getGmat(Jac,dXidu,fe.G);
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.basis(geoBasis);
|
||||
X.assign(Xnod * fe.basis(geoBasis));
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
@ -704,6 +705,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
|
||||
|
||||
std::vector<Matrix> dNxdu(m_basis.size());
|
||||
Matrix Xnod, Jac;
|
||||
double param[3] = { fe.u, fe.v, fe.w };
|
||||
Vec4 X;
|
||||
Vec3 normal;
|
||||
for (size_t l = 0; l < threadGrp[g][t].size() && ok; ++l)
|
||||
@ -776,17 +778,17 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
|
||||
if (gpar[0].size() > 1)
|
||||
{
|
||||
fe.xi = xg[k1];
|
||||
fe.u = gpar[0](k1+1,i1-p1+1);
|
||||
fe.u = param[0] = gpar[0](k1+1,i1-p1+1);
|
||||
}
|
||||
if (gpar[1].size() > 1)
|
||||
{
|
||||
fe.eta = xg[k2];
|
||||
fe.v = gpar[1](k2+1,i2-p2+1);
|
||||
fe.v = param[1] = gpar[1](k2+1,i2-p2+1);
|
||||
}
|
||||
if (gpar[2].size() > 1)
|
||||
{
|
||||
fe.zeta = xg[k3];
|
||||
fe.w = gpar[2](k3+1,i3-p3+1);
|
||||
fe.w = param[2] = gpar[2](k3+1,i3-p3+1);
|
||||
}
|
||||
|
||||
// Fetch basis function derivatives at current integration point
|
||||
|
||||
@ -34,20 +34,20 @@
|
||||
|
||||
|
||||
Fields* Fields::create (const ASMbase* pch, const RealArray& v,
|
||||
char basis, const char* name)
|
||||
char basis, int nf, const char* name)
|
||||
{
|
||||
#ifdef HAS_LRSPLINE
|
||||
const ASMu2Dmx* pu2mx = dynamic_cast<const ASMu2Dmx*>(pch);
|
||||
if (basis > 10 && pu2mx) return new LRSplineFields2Dmx(pu2mx,v,basis,name);
|
||||
|
||||
const ASMu2D* pu2 = dynamic_cast<const ASMu2D*>(pch);
|
||||
if (pu2) return new LRSplineFields2D(pu2,v,basis,name);
|
||||
if (pu2) return new LRSplineFields2D(pu2,v,basis,nf,name);
|
||||
|
||||
const ASMu3Dmx* pu3mx = dynamic_cast<const ASMu3Dmx*>(pch);
|
||||
if (basis > 10 && pu3mx) return new LRSplineFields3Dmx(pu3mx,v,basis,name);
|
||||
|
||||
const ASMu3D* pu3 = dynamic_cast<const ASMu3D*>(pch);
|
||||
if (pu3) return new LRSplineFields3D(pu3,v,basis,name);
|
||||
if (pu3) return new LRSplineFields3D(pu3,v,basis,nf,name);
|
||||
#endif
|
||||
const ASMs2DLag* pl2 = dynamic_cast<const ASMs2DLag*>(pch);
|
||||
if (pl2) return new LagrangeFields2D(pl2,v,basis,name);
|
||||
@ -56,7 +56,7 @@ Fields* Fields::create (const ASMbase* pch, const RealArray& v,
|
||||
if (basis > 10 && ps2mx) return new SplineFields2Dmx(ps2mx,v,basis,name);
|
||||
|
||||
const ASMs2D* ps2 = dynamic_cast<const ASMs2D*>(pch);
|
||||
if (ps2) return new SplineFields2D(ps2,v,basis,name);
|
||||
if (ps2) return new SplineFields2D(ps2,v,basis,nf,name);
|
||||
|
||||
const ASMs3DLag* pl3 = dynamic_cast<const ASMs3DLag*>(pch);
|
||||
if (pl3) return new LagrangeFields3D(pl3,v,basis,name);
|
||||
@ -65,7 +65,7 @@ Fields* Fields::create (const ASMbase* pch, const RealArray& v,
|
||||
if (basis > 10 && ps3mx) return new SplineFields3Dmx(ps3mx,v,basis,name);
|
||||
|
||||
const ASMs3D* ps3 = dynamic_cast<const ASMs3D*>(pch);
|
||||
if (ps3) return new SplineFields3D(ps3,v,basis,name);
|
||||
if (ps3) return new SplineFields3D(ps3,v,basis,nf,name);
|
||||
|
||||
return nullptr;
|
||||
}
|
||||
|
||||
@ -59,9 +59,10 @@ public:
|
||||
//! \param[in] pch The spline patch on which the field is to be defined on
|
||||
//! \param[in] v Array of nodal/control point field values
|
||||
//! \param[in] basis Basis to use from patch
|
||||
//! \param[in] nf Number of components in field
|
||||
//! \param[in] name Name of field
|
||||
static Fields* create(const ASMbase* pch, const RealArray& v,
|
||||
char basis = 1, const char* name = nullptr);
|
||||
char basis = 1, int nf = 0, const char* name = nullptr);
|
||||
|
||||
// Methods to evaluate the field
|
||||
//==============================
|
||||
|
||||
@ -965,7 +965,8 @@ bool ASMu2D::integrate (Integrand& integrand,
|
||||
fe.iel = MLGE[iel-1];
|
||||
Matrix dNdu, Xnod, Jac;
|
||||
Matrix3D d2Ndu2, Hess;
|
||||
Vec4 X;
|
||||
double param[3] = { 0.0, 0.0, 0.0 };
|
||||
Vec4 X(param);
|
||||
|
||||
// Get element area in the parameter space
|
||||
double dA = this->getParametricArea(iel);
|
||||
@ -1025,8 +1026,8 @@ bool ASMu2D::integrate (Integrand& integrand,
|
||||
fe.eta = xr[j];
|
||||
|
||||
// Parameter values of current integration point
|
||||
fe.u = redpar[0][i];
|
||||
fe.v = redpar[1][j];
|
||||
fe.u = param[0] = redpar[0][i];
|
||||
fe.v = param[1] = redpar[1][j];
|
||||
|
||||
// Compute basis function derivatives at current point
|
||||
Go::BasisDerivsSf spline;
|
||||
@ -1037,7 +1038,7 @@ bool ASMu2D::integrate (Integrand& integrand,
|
||||
fe.detJxW = utl::Jacobian(Jac,fe.dNdX,Xnod,dNdu);
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.N;
|
||||
X.assign(Xnod * fe.N);
|
||||
X.t = time.t;
|
||||
|
||||
// Compute the reduced integration terms of the integrand
|
||||
@ -1064,8 +1065,8 @@ bool ASMu2D::integrate (Integrand& integrand,
|
||||
fe.eta = xg[j];
|
||||
|
||||
// Parameter values of current integration point
|
||||
fe.u = gpar[0][i];
|
||||
fe.v = gpar[1][j];
|
||||
fe.u = param[0] = gpar[0][i];
|
||||
fe.v = param[1] = gpar[1][j];
|
||||
|
||||
// Compute basis function derivatives at current integration point
|
||||
if (integrand.getIntegrandType() & Integrand::SECOND_DERIVATIVES) {
|
||||
@ -1111,7 +1112,7 @@ bool ASMu2D::integrate (Integrand& integrand,
|
||||
#endif
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.N;
|
||||
X.assign(Xnod * fe.N);
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
@ -1271,6 +1272,7 @@ bool ASMu2D::integrate (Integrand& integrand,
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.N;
|
||||
X.u = elmPts[ip].data();
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
@ -1345,7 +1347,8 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex,
|
||||
size_t firstp = iit == firstBp.end() ? 0 : iit->second;
|
||||
|
||||
Matrix dNdu, Xnod, Jac;
|
||||
Vec4 X;
|
||||
double param[3] = { 0.0, 0.0, 0.0 };
|
||||
Vec4 X(param);
|
||||
Vec3 normal;
|
||||
|
||||
|
||||
@ -1401,8 +1404,8 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex,
|
||||
// of current integration point
|
||||
fe.xi = xg[i];
|
||||
fe.eta = xg[i];
|
||||
fe.u = gpar[0][i];
|
||||
fe.v = gpar[1][i];
|
||||
fe.u = param[0] = gpar[0][i];
|
||||
fe.v = param[1] = gpar[1][i];
|
||||
|
||||
// Evaluate basis function derivatives at current integration points
|
||||
Go::BasisDerivsSf spline;
|
||||
@ -1418,7 +1421,7 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex,
|
||||
if (edgeDir < 0) normal *= -1.0;
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.N;
|
||||
X.assign(Xnod * fe.N);
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
@ -1586,6 +1589,7 @@ bool ASMu2D::tesselate (ElementBlock& grid, const int* npe) const
|
||||
double v = vmin + (vmax-vmin)/(npe[1]-1)*iv;
|
||||
Go::Point pt;
|
||||
lrspline->point(pt, u,v, iel, iu!=npe[0]-1, iv!=npe[1]-1);
|
||||
grid.setParams(inod, u, v);
|
||||
for(int dim=0; dim<nsd; dim++)
|
||||
grid.setCoor(inod, dim, pt[dim]);
|
||||
inod++;
|
||||
|
||||
@ -76,7 +76,8 @@ LR::LRSplineSurface* ASMu2Dmx::getBasis (int basis)
|
||||
|
||||
bool ASMu2Dmx::write (std::ostream& os, int basis) const
|
||||
{
|
||||
return false;
|
||||
os << *m_basis[basis-1];
|
||||
return os.good();
|
||||
}
|
||||
|
||||
|
||||
@ -288,7 +289,8 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
|
||||
fe.iel = MLGE[geoEl-1];
|
||||
std::vector<Matrix> dNxdu(m_basis.size());
|
||||
Matrix Xnod, Jac;
|
||||
Vec4 X;
|
||||
double param[3] = { 0.0, 0.0, 0.0 };
|
||||
Vec4 X(param);
|
||||
std::vector<Matrix3D> d2Nxdu2(m_basis.size());
|
||||
Matrix3D Hess;
|
||||
|
||||
@ -337,8 +339,8 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
|
||||
fe.eta = xg[j];
|
||||
|
||||
// Parameter values of current integration point
|
||||
fe.u = gpar[0][i];
|
||||
fe.v = gpar[1][j];
|
||||
fe.u = param[0] = gpar[0][i];
|
||||
fe.v = param[1] = gpar[1][j];
|
||||
|
||||
// Compute basis function derivatives at current integration point
|
||||
if (use2ndDer)
|
||||
@ -375,7 +377,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand,
|
||||
}
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.basis(geoBasis);
|
||||
X.assign(Xnod * fe.basis(geoBasis));
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
@ -451,7 +453,8 @@ bool ASMu2Dmx::integrate (Integrand& integrand, int lIndex,
|
||||
|
||||
std::vector<Matrix> dNxdu(m_basis.size());
|
||||
Matrix Xnod, Jac;
|
||||
Vec4 X;
|
||||
double param[3] = { 0.0, 0.0, 0.0 };
|
||||
Vec4 X(param);
|
||||
Vec3 normal;
|
||||
|
||||
// === Assembly loop over all elements on the patch edge =====================
|
||||
@ -516,8 +519,8 @@ bool ASMu2Dmx::integrate (Integrand& integrand, int lIndex,
|
||||
// of current integration point
|
||||
fe.xi = xg[i];
|
||||
fe.eta = xg[i];
|
||||
fe.u = gpar[0][i];
|
||||
fe.v = gpar[1][i];
|
||||
fe.u = param[0] = gpar[0][i];
|
||||
fe.v = param[1] = gpar[1][i];
|
||||
|
||||
// Evaluate basis function derivatives at current integration points
|
||||
std::vector<Go::BasisDerivsSf> splinex(m_basis.size());
|
||||
@ -538,7 +541,7 @@ bool ASMu2Dmx::integrate (Integrand& integrand, int lIndex,
|
||||
normal *= -1.0;
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.basis(geoBasis);
|
||||
X.assign(Xnod * fe.basis(geoBasis));
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
|
||||
@ -1018,7 +1018,8 @@ bool ASMu3D::integrate (Integrand& integrand,
|
||||
Matrix dNdu, Xnod, Jac;
|
||||
Matrix3D d2Ndu2, Hess;
|
||||
double dXidu[3];
|
||||
Vec4 X;
|
||||
double param[3] = { 0.0, 0.0, 0.0 };
|
||||
Vec4 X(param);
|
||||
// Get element volume in the parameter space
|
||||
double du = el->umax() - el->umin();
|
||||
double dv = el->vmax() - el->vmin();
|
||||
@ -1163,9 +1164,9 @@ bool ASMu3D::integrate (Integrand& integrand,
|
||||
fe.zeta = xg[k];
|
||||
|
||||
// Parameter values of current integration point
|
||||
fe.u = gpar[0][i];
|
||||
fe.v = gpar[1][j];
|
||||
fe.w = gpar[2][k];
|
||||
fe.u = param[0] = gpar[0][i];
|
||||
fe.v = param[1] = gpar[1][j];
|
||||
fe.w = param[2] = gpar[2][k];
|
||||
|
||||
// Extract bezier basis functions
|
||||
B.fillColumn(1, BN.getColumn(ig));
|
||||
@ -1240,7 +1241,7 @@ bool ASMu3D::integrate (Integrand& integrand,
|
||||
utl::getGmat(Jac,dXidu,fe.G);
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.N;
|
||||
X.assign(Xnod * fe.N);
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
@ -1346,7 +1347,8 @@ bool ASMu3D::integrate (Integrand& integrand, int lIndex,
|
||||
fe.w = gpar[2](1);
|
||||
|
||||
Matrix dNdu, Xnod, Jac;
|
||||
Vec4 X;
|
||||
double param[3] = { fe.u, fe.v, fe.w };
|
||||
Vec4 X(param);
|
||||
Vec3 normal;
|
||||
double dXidu[3];
|
||||
|
||||
@ -1404,17 +1406,17 @@ bool ASMu3D::integrate (Integrand& integrand, int lIndex,
|
||||
if (gpar[0].size() > 1)
|
||||
{
|
||||
fe.xi = xg[k1];
|
||||
fe.u = gpar[0](k1+1);
|
||||
fe.u = param[0] = gpar[0](k1+1);
|
||||
}
|
||||
if (gpar[1].size() > 1)
|
||||
{
|
||||
fe.eta = xg[k2];
|
||||
fe.v = gpar[1](k2+1);
|
||||
fe.v = param[1] = gpar[1](k2+1);
|
||||
}
|
||||
if (gpar[2].size() > 1)
|
||||
{
|
||||
fe.zeta = xg[k3];
|
||||
fe.w = gpar[2](k3+1);
|
||||
fe.w = param[2] = gpar[2](k3+1);
|
||||
}
|
||||
|
||||
// Fetch basis function derivatives at current integration point
|
||||
@ -1433,7 +1435,7 @@ bool ASMu3D::integrate (Integrand& integrand, int lIndex,
|
||||
utl::getGmat(Jac,dXidu,fe.G);
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.N;
|
||||
X.assign(Xnod * fe.N);
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
@ -1753,6 +1755,7 @@ bool ASMu3D::tesselate (ElementBlock& grid, const int* npe) const
|
||||
double w = wmin + (wmax-wmin)/(npe[2]-1)*iw;
|
||||
Go::Point pt;
|
||||
lrspline->point(pt, u,v,w, iel, iu!=npe[0]-1, iv!=npe[1]-1, iw!=npe[2]-1);
|
||||
grid.setParams(inod, u, v, w);
|
||||
for(int dim=0; dim<nsd; dim++)
|
||||
grid.setCoor(inod, dim, pt[dim]);
|
||||
inod++;
|
||||
|
||||
@ -358,7 +358,8 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
|
||||
std::vector<Matrix3D> d2Nxdu2(m_basis.size());
|
||||
Matrix3D Hess;
|
||||
double dXidu[3];
|
||||
Vec4 X;
|
||||
double param[3] = { 0.0, 0.0, 0.0 };
|
||||
Vec4 X(param);
|
||||
// Get element volume in the parameter space
|
||||
double du = el->umax() - el->umin();
|
||||
double dv = el->vmax() - el->vmin();
|
||||
@ -468,9 +469,9 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
|
||||
fe.zeta = xg[k];
|
||||
|
||||
// Parameter values of current integration point
|
||||
fe.u = gpar[0][i];
|
||||
fe.v = gpar[1][j];
|
||||
fe.w = gpar[2][k];
|
||||
fe.u = param[0] = gpar[0][i];
|
||||
fe.v = param[1] = gpar[1][j];
|
||||
fe.w = param[2] = gpar[2][k];
|
||||
|
||||
// Extract bezier basis functions
|
||||
for (size_t b = 0; b < m_basis.size(); ++b) {
|
||||
@ -512,7 +513,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand,
|
||||
utl::getGmat(Jac,dXidu,fe.G);
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.basis(geoBasis);
|
||||
X.assign(Xnod * fe.basis(geoBasis));
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
@ -613,7 +614,9 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex,
|
||||
|
||||
std::vector<Matrix> dNxdu(m_basis.size());
|
||||
Matrix Xnod, Jac;
|
||||
Vec4 X;
|
||||
|
||||
double param[3] = { fe.u, fe.v, fe.w };
|
||||
Vec4 X(param);
|
||||
Vec3 normal;
|
||||
double dXidu[3];
|
||||
|
||||
@ -671,17 +674,17 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex,
|
||||
if (gpar[0].size() > 1)
|
||||
{
|
||||
fe.xi = xg[k1];
|
||||
fe.u = gpar[0](k1+1);
|
||||
fe.u = param[0] = gpar[0](k1+1);
|
||||
}
|
||||
if (gpar[1].size() > 1)
|
||||
{
|
||||
fe.eta = xg[k2];
|
||||
fe.v = gpar[1](k2+1);
|
||||
fe.v = param[1] = gpar[1](k2+1);
|
||||
}
|
||||
if (gpar[2].size() > 1)
|
||||
{
|
||||
fe.zeta = xg[k3];
|
||||
fe.w = gpar[2](k3+1);
|
||||
fe.w = param[2] = gpar[2](k3+1);
|
||||
}
|
||||
|
||||
// Fetch basis function derivatives at current integration point
|
||||
@ -703,7 +706,7 @@ bool ASMu3Dmx::integrate (Integrand& integrand, int lIndex,
|
||||
utl::getGmat(Jac,dXidu,fe.G);
|
||||
|
||||
// Cartesian coordinates of current integration point
|
||||
X = Xnod * fe.basis(geoBasis);
|
||||
X.assign(Xnod * fe.basis(geoBasis));
|
||||
X.t = time.t;
|
||||
|
||||
// Evaluate the integrand and accumulate element contributions
|
||||
|
||||
@ -73,6 +73,14 @@ double LRSplineField2D::valueFE (const FiniteElement& fe) const
|
||||
|
||||
double LRSplineField2D::valueCoor (const Vec4& x) const
|
||||
{
|
||||
if (x.u) {
|
||||
FiniteElement fe;
|
||||
fe.u = x.u[0];
|
||||
fe.v = x.u[1];
|
||||
|
||||
return this->valueFE(fe);
|
||||
}
|
||||
|
||||
assert(0);
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
@ -73,6 +73,15 @@ double LRSplineField3D::valueFE (const FiniteElement& fe) const
|
||||
|
||||
double LRSplineField3D::valueCoor (const Vec4& x) const
|
||||
{
|
||||
if (x.u) {
|
||||
FiniteElement fe;
|
||||
fe.u = x.u[0];
|
||||
fe.v = x.u[1];
|
||||
fe.w = x.u[2];
|
||||
|
||||
return this->valueFE(fe);
|
||||
}
|
||||
|
||||
assert(0);
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
@ -22,7 +22,7 @@
|
||||
|
||||
LRSplineFields2D::LRSplineFields2D (const ASMu2D* patch,
|
||||
const RealArray& v, char nbasis,
|
||||
const char* name)
|
||||
int nnf, const char* name)
|
||||
: Fields(name), basis(patch->getBasis(nbasis)), surf(patch->getSurface())
|
||||
{
|
||||
nno = basis->nBasisFunctions();
|
||||
@ -33,17 +33,19 @@ LRSplineFields2D::LRSplineFields2D (const ASMu2D* patch,
|
||||
ofs += patch->getNoNodes(i)*patch->getNoFields(i);
|
||||
auto vit = v.begin()+ofs;
|
||||
|
||||
nf = 2;
|
||||
if (nnf == 0)
|
||||
nnf = 2;
|
||||
nf = nnf;
|
||||
int nfc = patch->getNoFields(nbasis);
|
||||
values.resize(nno*nf);
|
||||
int ndof = nfc*nno;
|
||||
auto end = v.size() > ofs+ndof ? vit+ndof : v.end();
|
||||
if (nfc == 2)
|
||||
if (nfc == nf)
|
||||
std::copy(vit,end,values.begin());
|
||||
else
|
||||
for (size_t i = 0; i < nno && vit != end; ++i, vit += nfc)
|
||||
for (size_t j = 0; j < 2; ++j)
|
||||
values[2*i+j] = *(vit+j);
|
||||
for (size_t j = 0; j < nf; ++j)
|
||||
values[nf*i+j] = *(vit+j);
|
||||
}
|
||||
|
||||
|
||||
@ -69,11 +71,11 @@ bool LRSplineFields2D::valueFE (const FiniteElement& fe, Vector& vals) const
|
||||
basis->computeBasis(fe.u,fe.v,spline,iel);
|
||||
|
||||
// Evaluate the solution field at the given point
|
||||
Matrix Vnod(2, elm->nBasisFunctions());
|
||||
Matrix Vnod(nf, elm->nBasisFunctions());
|
||||
size_t i = 1;
|
||||
for (auto it = elm->constSupportBegin();
|
||||
it != elm->constSupportEnd(); ++it, ++i)
|
||||
for (size_t j = 1; j <= 2; ++j)
|
||||
for (size_t j = 1; j <= nf; ++j)
|
||||
Vnod(j,i) = values((*it)->getId()*2+j);
|
||||
|
||||
Vnod.multiply(spline.basisValues,vals); // vals = Vnod * basisValues
|
||||
@ -84,6 +86,14 @@ bool LRSplineFields2D::valueFE (const FiniteElement& fe, Vector& vals) const
|
||||
|
||||
bool LRSplineFields2D::valueCoor (const Vec4& x, Vector& vals) const
|
||||
{
|
||||
if (x.u) {
|
||||
FiniteElement fe;
|
||||
fe.u = x.u[0];
|
||||
fe.v = x.u[1];
|
||||
|
||||
return this->valueFE(fe, vals);
|
||||
}
|
||||
|
||||
assert(0);
|
||||
return false;
|
||||
}
|
||||
@ -137,20 +147,20 @@ bool LRSplineFields2D::gradFE (const FiniteElement& fe, Matrix& grad) const
|
||||
}
|
||||
dNdX.multiply(dNdu,Jac); // dNdX = dNdu * Jac
|
||||
|
||||
Vnod.resize(2, nbf);
|
||||
Vnod.resize(nf, nbf);
|
||||
i = 1;
|
||||
for (auto it = belm->constSupportBegin();
|
||||
it != belm->constSupportEnd(); ++it, ++i)
|
||||
for (size_t j = 1; j <= 2; ++j)
|
||||
Vnod(j,i) = values((*it)->getId()*2+j);
|
||||
for (size_t j = 1; j <= nf; ++j)
|
||||
Vnod(j,i) = values((*it)->getId()*nf+j);
|
||||
}
|
||||
else {
|
||||
Vnod.resize(2, nen);
|
||||
Vnod.resize(nf, nen);
|
||||
i = 1;
|
||||
for (auto it = elm->constSupportBegin();
|
||||
it != elm->constSupportEnd(); ++it, ++i)
|
||||
for (size_t j = 1; j <= 2; ++j)
|
||||
Vnod(j,i) = values((*it)->getId()*2+j);
|
||||
for (size_t j = 1; j <= nf; ++j)
|
||||
Vnod(j,i) = values((*it)->getId()*nf+j);
|
||||
}
|
||||
|
||||
grad.multiply(Vnod,dNdX); // grad = Vnod * dNdX
|
||||
@ -193,12 +203,12 @@ bool LRSplineFields2D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
|
||||
d2Ndu2(n,2,2) = spline2.basisDerivs_vv[n-1];
|
||||
}
|
||||
|
||||
Vnod.resize(2, nen);
|
||||
Vnod.resize(nf, nen);
|
||||
i = 1;
|
||||
for (auto it = elm->constSupportBegin();
|
||||
it != elm->constSupportEnd(); ++it, ++i)
|
||||
for (size_t j = 1; j <= 2; ++j)
|
||||
Vnod(j,i) = values((*it)->getId()*2+j);
|
||||
for (size_t j = 1; j <= nf; ++j)
|
||||
Vnod(j,i) = values((*it)->getId()*nf+j);
|
||||
}
|
||||
else {
|
||||
surf->computeBasis(fe.u,fe.v,spline,iel);
|
||||
@ -231,12 +241,12 @@ bool LRSplineFields2D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
|
||||
d2Ndu2(n,2,2) = spline2.basisDerivs_vv[n-1];
|
||||
}
|
||||
|
||||
Vnod.resize(2, nbf);
|
||||
Vnod.resize(nf, nbf);
|
||||
i = 1;
|
||||
for (auto it = belm->constSupportBegin();
|
||||
it != belm->constSupportEnd(); ++it, ++i)
|
||||
for (size_t j = 1; j <= 2; ++j)
|
||||
Vnod(j,i) = values((*it)->getId()*2+j);
|
||||
for (size_t j = 1; j <= nf; ++j)
|
||||
Vnod(j,i) = values((*it)->getId()*nf+j);
|
||||
}
|
||||
|
||||
return H.multiply(Vnod,d2Ndu2);
|
||||
|
||||
@ -37,9 +37,11 @@ public:
|
||||
//! \param[in] patch The spline patch on which the field is to be defined
|
||||
//! \param[in] v Array of control point field values
|
||||
//! \param[in] basis Basis to use from patch
|
||||
//! \param[in] nf Number of components to for field
|
||||
//! \param[in] name Name of spline field
|
||||
LRSplineFields2D(const ASMu2D* patch, const RealArray& v,
|
||||
char basis = 1, const char* name = nullptr);
|
||||
char basis = 1, int nf = 0, const char* name = nullptr);
|
||||
|
||||
//! \brief Empty destructor.
|
||||
virtual ~LRSplineFields2D() {}
|
||||
|
||||
|
||||
@ -22,7 +22,7 @@
|
||||
|
||||
LRSplineFields3D::LRSplineFields3D (const ASMu3D* patch,
|
||||
const RealArray& v, char nbasis,
|
||||
const char* name)
|
||||
int nnf, const char* name)
|
||||
: Fields(name), basis(patch->getBasis(nbasis)), vol(patch->getVolume())
|
||||
{
|
||||
nno = basis->nBasisFunctions();
|
||||
@ -33,17 +33,19 @@ LRSplineFields3D::LRSplineFields3D (const ASMu3D* patch,
|
||||
ofs += patch->getNoNodes(i)*patch->getNoFields(i);
|
||||
auto vit = v.begin()+ofs;
|
||||
|
||||
nf = 3;
|
||||
if (nnf == 0)
|
||||
nnf = 3;
|
||||
nf = nnf;
|
||||
int nfc = patch->getNoFields(nbasis);
|
||||
values.resize(nno*nf);
|
||||
int ndof = nfc*nno;
|
||||
auto end = v.size() > ofs+ndof ? vit+ndof : v.end();
|
||||
if (nfc == 3)
|
||||
if (nfc == nf)
|
||||
std::copy(vit,end,values.begin());
|
||||
else
|
||||
for (size_t i = 0; i < nno && vit != end; ++i, vit += nfc)
|
||||
for (size_t j = 0; j < 3; ++j)
|
||||
values[3*i+j] = *(vit+j);
|
||||
for (size_t j = 0; j < nf; ++j)
|
||||
values[nf*i+j] = *(vit+j);
|
||||
}
|
||||
|
||||
|
||||
@ -69,12 +71,12 @@ bool LRSplineFields3D::valueFE (const FiniteElement& fe, Vector& vals) const
|
||||
basis->computeBasis(fe.u,fe.v,fe.w,spline,iel);
|
||||
|
||||
// Evaluate the solution field at the given point
|
||||
Matrix Vnod(3, elm->nBasisFunctions());
|
||||
Matrix Vnod(nf, elm->nBasisFunctions());
|
||||
size_t i = 1;
|
||||
for (auto it = elm->constSupportBegin();
|
||||
it != elm->constSupportEnd(); ++it, ++i)
|
||||
for (size_t j = 1; j <= 3; ++j)
|
||||
Vnod(j,i) = values((*it)->getId()*3+j);
|
||||
for (size_t j = 1; j <= nf; ++j)
|
||||
Vnod(j,i) = values((*it)->getId()*nf+j);
|
||||
|
||||
Vnod.multiply(spline.basisValues,vals); // vals = Vnod * basisValues
|
||||
|
||||
@ -84,6 +86,15 @@ bool LRSplineFields3D::valueFE (const FiniteElement& fe, Vector& vals) const
|
||||
|
||||
bool LRSplineFields3D::valueCoor (const Vec4& x, Vector& vals) const
|
||||
{
|
||||
if (x.u) {
|
||||
FiniteElement fe;
|
||||
fe.u = x.u[0];
|
||||
fe.v = x.u[1];
|
||||
fe.w = x.u[2];
|
||||
|
||||
return this->valueFE(fe, vals);
|
||||
}
|
||||
|
||||
assert(0);
|
||||
return false;
|
||||
}
|
||||
@ -139,20 +150,20 @@ bool LRSplineFields3D::gradFE (const FiniteElement& fe, Matrix& grad) const
|
||||
}
|
||||
dNdX.multiply(dNdu,Jac); // dNdX = dNdu * Jac
|
||||
|
||||
Vnod.resize(3, nbf);
|
||||
Vnod.resize(nf, nbf);
|
||||
i = 1;
|
||||
for (auto it = belm->constSupportBegin();
|
||||
it != belm->constSupportEnd(); ++it, ++i)
|
||||
for (size_t j = 1; j <= 3; ++j)
|
||||
Vnod(j,i) = values((*it)->getId()*3+j);
|
||||
for (size_t j = 1; j <= nf; ++j)
|
||||
Vnod(j,i) = values((*it)->getId()*nf+j);
|
||||
}
|
||||
else {
|
||||
Vnod.resize(3, nen);
|
||||
Vnod.resize(nf, nen);
|
||||
i = 1;
|
||||
for (auto it = elm->constSupportBegin();
|
||||
it != elm->constSupportEnd(); ++it, ++i)
|
||||
for (size_t j = 1; j <= 3; ++j)
|
||||
Vnod(j,i) = values((*it)->getId()*3+j);
|
||||
for (size_t j = 1; j <= nf; ++j)
|
||||
Vnod(j,i) = values((*it)->getId()*nf+j);
|
||||
}
|
||||
|
||||
grad.multiply(Vnod,dNdX); // grad = Vnod * dNdX
|
||||
@ -199,12 +210,12 @@ bool LRSplineFields3D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
|
||||
d2Ndu2(n,3,3) = spline2.basisDerivs_ww[n-1];
|
||||
}
|
||||
|
||||
Vnod.resize(3, nen);
|
||||
Vnod.resize(nf, nen);
|
||||
i = 1;
|
||||
for (auto it = elm->constSupportBegin();
|
||||
it != elm->constSupportEnd(); ++it, ++i)
|
||||
for (size_t j = 1; j <= 3; ++j)
|
||||
Vnod(j,i) = values((*it)->getId()*3+j);
|
||||
for (size_t j = 1; j <= nf; ++j)
|
||||
Vnod(j,i) = values((*it)->getId()*nf+j);
|
||||
}
|
||||
else {
|
||||
vol->computeBasis(fe.u,fe.v,fe.w,spline,iel);
|
||||
@ -242,12 +253,12 @@ bool LRSplineFields3D::hessianFE (const FiniteElement& fe, Matrix3D& H) const
|
||||
d2Ndu2(n,3,3) = spline2.basisDerivs_ww[n-1];
|
||||
}
|
||||
|
||||
Vnod.resize(3, nbf);
|
||||
Vnod.resize(nf, nbf);
|
||||
i = 1;
|
||||
for (auto it = belm->constSupportBegin();
|
||||
it != belm->constSupportEnd(); ++it, ++i)
|
||||
for (size_t j = 1; j <= 3; ++j)
|
||||
Vnod(j,i) = values((*it)->getId()*3+j);
|
||||
for (size_t j = 1; j <= nf; ++j)
|
||||
Vnod(j,i) = values((*it)->getId()*nf+j);
|
||||
}
|
||||
|
||||
return H.multiply(Vnod,d2Ndu2);
|
||||
|
||||
@ -37,9 +37,11 @@ public:
|
||||
//! \param[in] patch The spline patch on which the field is to be defined
|
||||
//! \param[in] v Array of control point field values
|
||||
//! \param[in] basis Basis to use from patch
|
||||
//! \param[in] nf Number of components to for field
|
||||
//! \param[in] name Name of spline field
|
||||
LRSplineFields3D(const ASMu3D* patch, const RealArray& v,
|
||||
char basis = 1, const char* name = nullptr);
|
||||
char basis = 1, int nf = 0, const char* name = nullptr);
|
||||
|
||||
//! \brief Empty destructor.
|
||||
virtual ~LRSplineFields3D() {}
|
||||
|
||||
|
||||
@ -81,20 +81,25 @@ double SplineField2D::valueFE (const FiniteElement& fe) const
|
||||
|
||||
double SplineField2D::valueCoor (const Vec4& x) const
|
||||
{
|
||||
// use with caution
|
||||
Go::Point pt(3), clopt(3);
|
||||
pt[0] = x[0];
|
||||
pt[1] = x[1];
|
||||
pt[2] = x[2];
|
||||
double clo_u, clo_v, dist;
|
||||
#pragma omp critical
|
||||
surf->closestPoint(pt, clo_u, clo_v, clopt, dist, 1e-5);
|
||||
|
||||
FiniteElement fe;
|
||||
fe.u = clopt[0];
|
||||
fe.v = clopt[1];
|
||||
if (x.u) {
|
||||
fe.u = x.u[0];
|
||||
fe.v = x.u[1];
|
||||
}
|
||||
else {
|
||||
// use with caution
|
||||
Go::Point pt(3), clopt(3);
|
||||
pt[0] = x[0];
|
||||
pt[1] = x[1];
|
||||
pt[2] = x[2];
|
||||
double clo_u, clo_v, dist;
|
||||
#pragma omp critical
|
||||
surf->closestPoint(pt, clo_u, clo_v, clopt, dist, 1e-10);
|
||||
fe.u = clo_u;
|
||||
fe.v = clo_v;
|
||||
}
|
||||
|
||||
return valueFE(fe);
|
||||
return this->valueFE(fe);
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -82,8 +82,28 @@ double SplineField3D::valueFE (const FiniteElement& fe) const
|
||||
|
||||
double SplineField3D::valueCoor (const Vec4& x) const
|
||||
{
|
||||
// Not implemented yet
|
||||
return 0.0;
|
||||
FiniteElement fe;
|
||||
if (x.u) {
|
||||
fe.u = x.u[0];
|
||||
fe.v = x.u[1];
|
||||
fe.w = x.u[2];
|
||||
}
|
||||
else {
|
||||
// use with caution
|
||||
Go::Point pt(3), clopt(3);
|
||||
pt[0] = x[0];
|
||||
pt[1] = x[1];
|
||||
pt[2] = x[2];
|
||||
double clo_u, clo_v, clo_w, dist;
|
||||
#pragma omp critical
|
||||
vol->closestPoint(pt, clo_u, clo_v, clo_w, clopt, dist, 1e-5);
|
||||
|
||||
fe.u = clo_u;
|
||||
fe.v = clo_v;
|
||||
fe.w = clo_w;
|
||||
}
|
||||
|
||||
return this->valueFE(fe);
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -23,7 +23,7 @@
|
||||
|
||||
SplineFields2D::SplineFields2D (const ASMs2D* patch,
|
||||
const RealArray& v, char nbasis,
|
||||
const char* name)
|
||||
int nnf, const char* name)
|
||||
: Fields(name), basis(patch->getBasis(nbasis)), surf(patch->getSurface())
|
||||
{
|
||||
const int n1 = basis->numCoefs_u();
|
||||
@ -40,17 +40,20 @@ SplineFields2D::SplineFields2D (const ASMs2D* patch,
|
||||
auto vit = v.begin()+ofs;
|
||||
|
||||
// Ensure the values array has compatible length, pad with zeros if necessary
|
||||
nf = 2;
|
||||
if (nnf == 0)
|
||||
nnf = 2;
|
||||
|
||||
nf = nnf;
|
||||
int nfc = patch->getNoFields(nbasis);
|
||||
values.resize(nno*nf);
|
||||
int ndof = nfc*nno;
|
||||
auto end = v.size() > ofs+ndof ? vit+ndof : v.end();
|
||||
if (nfc == 2)
|
||||
if (nfc == nf)
|
||||
std::copy(vit,end,values.begin());
|
||||
else
|
||||
for (size_t i = 0; i < nno && vit != end; ++i, vit += nfc)
|
||||
for (size_t j = 0; j < 2; ++j)
|
||||
values[2*i+j] = *(vit+j);
|
||||
for (size_t j = 0; j < nf; ++j)
|
||||
values[nf*i+j] = *(vit+j);
|
||||
}
|
||||
|
||||
|
||||
@ -89,8 +92,26 @@ bool SplineFields2D::valueFE (const FiniteElement& fe, Vector& vals) const
|
||||
|
||||
bool SplineFields2D::valueCoor (const Vec4& x, Vector& vals) const
|
||||
{
|
||||
// Not implemented yet
|
||||
return false;
|
||||
FiniteElement fe;
|
||||
if (x.u) {
|
||||
fe.u = x.u[0];
|
||||
fe.v = x.u[1];
|
||||
}
|
||||
else {
|
||||
// use with caution
|
||||
Go::Point pt(3), clopt(3);
|
||||
pt[0] = x[0];
|
||||
pt[1] = x[1];
|
||||
pt[2] = x[2];
|
||||
double clo_u, clo_v, dist;
|
||||
#pragma omp critical
|
||||
surf->closestPoint(pt, clo_u, clo_v, clopt, dist, 1e-5);
|
||||
|
||||
fe.u = clo_u;
|
||||
fe.v = clo_v;
|
||||
}
|
||||
|
||||
return this->valueFE(fe, vals);
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -37,9 +37,20 @@ public:
|
||||
//! \param[in] patch The spline patch on which the field is to be defined
|
||||
//! \param[in] v Array of control point field values
|
||||
//! \param[in] basis Basis to use from patch
|
||||
//! \param[in] nf Number of components in field
|
||||
//! \param[in] name Name of spline field
|
||||
SplineFields2D(const ASMs2D* patch, const RealArray& v,
|
||||
char basis = 1, const char* name = nullptr);
|
||||
char basis = 1, int nf = 0, const char* name = nullptr);
|
||||
|
||||
//! \brief Construct directly from surface.
|
||||
//! \param[in] srf The spline surface to use
|
||||
//! \param[in] v Array of control point field values
|
||||
//! \param[in] ncmp Number of field components
|
||||
//! \param[in] name Name of spline field
|
||||
SplineFields2D(const Go::SplineSurface* srf, const RealArray& v, int ncmp,
|
||||
const char* name = nullptr);
|
||||
|
||||
|
||||
//! \brief Empty destructor.
|
||||
virtual ~SplineFields2D() {}
|
||||
|
||||
|
||||
@ -23,7 +23,7 @@
|
||||
|
||||
SplineFields3D::SplineFields3D (const ASMs3D* patch,
|
||||
const RealArray& v, char nbasis,
|
||||
const char* name)
|
||||
int nnf, const char* name)
|
||||
: Fields(name), basis(patch->getBasis(nbasis)), vol(patch->getVolume())
|
||||
{
|
||||
const int n1 = basis->numCoefs(0);
|
||||
@ -41,17 +41,20 @@ SplineFields3D::SplineFields3D (const ASMs3D* patch,
|
||||
ofs += patch->getNoNodes(i)*patch->getNoFields(i);
|
||||
auto vit = v.begin()+ofs;
|
||||
|
||||
nf = 3;
|
||||
if (nnf == 0)
|
||||
nnf = 3;
|
||||
|
||||
nf = nnf;
|
||||
int nfc = patch->getNoFields(nbasis);
|
||||
values.resize(nno*nf);
|
||||
int ndof = nfc*nno;
|
||||
auto end = v.size() > ofs+ndof ? vit+ndof : v.end();
|
||||
if (nfc == 3)
|
||||
if (nfc == nf)
|
||||
std::copy(vit,end,values.begin());
|
||||
else
|
||||
for (size_t i = 0; i < nno && vit != end; ++i, vit += nfc)
|
||||
for (size_t j = 0; j < 3; ++j)
|
||||
values[3*i+j] = *(vit+j);
|
||||
for (size_t j = 0; j < nf; ++j)
|
||||
values[nf*i+j] = *(vit+j);
|
||||
}
|
||||
|
||||
|
||||
@ -90,8 +93,28 @@ bool SplineFields3D::valueFE (const FiniteElement& fe, Vector& vals) const
|
||||
|
||||
bool SplineFields3D::valueCoor (const Vec4& x, Vector& vals) const
|
||||
{
|
||||
// Not implemented yet
|
||||
return false;
|
||||
FiniteElement fe;
|
||||
if (x.u) {
|
||||
fe.u = x.u[0];
|
||||
fe.v = x.u[1];
|
||||
fe.w = x.u[2];
|
||||
}
|
||||
else {
|
||||
// use with caution
|
||||
Go::Point pt(3), clopt(3);
|
||||
pt[0] = x[0];
|
||||
pt[1] = x[1];
|
||||
pt[2] = x[2];
|
||||
double clo_u, clo_v, clo_w, dist;
|
||||
#pragma omp critical
|
||||
vol->closestPoint(pt, clo_u, clo_v, clo_w, clopt, dist, 1e-5);
|
||||
|
||||
fe.u = clo_u;
|
||||
fe.v = clo_v;
|
||||
fe.w = clo_w;
|
||||
}
|
||||
|
||||
return this->valueFE(fe, vals);
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -37,9 +37,10 @@ public:
|
||||
//! \param[in] patch The spline patch on which the field is to be defined
|
||||
//! \param[in] v Array of control point field values
|
||||
//! \param[in] basis Basis to use from patch
|
||||
//! \param[in] nf Number of components in field
|
||||
//! \param[in] name Name of spline field
|
||||
SplineFields3D(const ASMs3D* patch, const RealArray& v,
|
||||
char basis = 1, const char* name = nullptr);
|
||||
char basis = 1, int nf = 0, const char* name = nullptr);
|
||||
//! \brief Empty destructor.
|
||||
virtual ~SplineFields3D() {}
|
||||
|
||||
|
||||
@ -810,6 +810,8 @@ bool SIMbase::assembleSystem (const TimeDomain& time, const Vectors& prevSol,
|
||||
<< std::endl;
|
||||
ok &= this->initBodyLoad(lp);
|
||||
ok &= this->extractPatchSolution(it->second,prevSol,lp-1);
|
||||
if (mySol)
|
||||
mySol->initPatch(pch->idx);
|
||||
ok &= pch->integrate(*it->second,sysQ,time);
|
||||
}
|
||||
else
|
||||
@ -826,6 +828,8 @@ bool SIMbase::assembleSystem (const TimeDomain& time, const Vectors& prevSol,
|
||||
<< std::endl;
|
||||
ok &= this->initBodyLoad(lp);
|
||||
ok &= this->extractPatchSolution(it->second,prevSol,k);
|
||||
if (mySol)
|
||||
mySol->initPatch(myModel[k]->idx);
|
||||
ok &= myModel[k]->integrate(*it->second,sysQ,time);
|
||||
}
|
||||
}
|
||||
@ -846,6 +850,9 @@ bool SIMbase::assembleSystem (const TimeDomain& time, const Vectors& prevSol,
|
||||
break;
|
||||
}
|
||||
|
||||
if (mySol)
|
||||
mySol->initPatch(pch->idx);
|
||||
|
||||
for (p2 = myProps.begin(); p2 != myProps.end() && ok; ++p2)
|
||||
if (p2->pcode == Property::MATERIAL && p->patch == p2->patch)
|
||||
if (!(ok = this->initMaterial(p2->pindx)))
|
||||
@ -1331,6 +1338,8 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
|
||||
norm->getProjection(k).clear();
|
||||
else
|
||||
this->extractPatchSolution(ssol[k],norm->getProjection(k),lp-1,nCmp,1);
|
||||
if (mySol)
|
||||
mySol->initPatch(pch->idx);
|
||||
ok &= pch->integrate(*norm,globalNorm,time);
|
||||
}
|
||||
else
|
||||
@ -1347,6 +1356,8 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
|
||||
norm->getProjection(k).clear();
|
||||
else
|
||||
this->extractPatchSolution(ssol[k],norm->getProjection(k),i,nCmp,1);
|
||||
if (mySol)
|
||||
mySol->initPatch(myModel[i]->idx);
|
||||
ok &= myModel[i]->integrate(*norm,globalNorm,time);
|
||||
lp = i+1;
|
||||
}
|
||||
@ -1362,6 +1373,8 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
|
||||
{
|
||||
if (p->patch != lp)
|
||||
ok = this->extractPatchSolution(psol,p->patch-1);
|
||||
if (mySol)
|
||||
mySol->initPatch(pch->idx);
|
||||
ok &= pch->integrate(*norm,p->lindx,globalNorm,time);
|
||||
lp = p->patch;
|
||||
}
|
||||
@ -1373,6 +1386,8 @@ bool SIMbase::solutionNorms (const TimeDomain& time,
|
||||
{
|
||||
if (p->patch != lp)
|
||||
ok = this->extractPatchSolution(psol,p->patch-1);
|
||||
if (mySol)
|
||||
mySol->initPatch(pch->idx);
|
||||
ok &= pch->integrateEdge(*norm,p->lindx,globalNorm,time);
|
||||
lp = p->patch;
|
||||
}
|
||||
|
||||
@ -547,6 +547,8 @@ int SIMoutput::writeGlvS1 (const Vector& psol, int iStep, int& nBlock,
|
||||
return -1;
|
||||
|
||||
myModel[i]->filterResults(field,myVtf->getBlock(++geomID));
|
||||
if (mySol)
|
||||
mySol->initPatch(myModel[i]->idx);
|
||||
|
||||
if (!scalarOnly && (nVcomp > 1 || !pvecName))
|
||||
{
|
||||
@ -576,7 +578,7 @@ int SIMoutput::writeGlvS1 (const Vector& psol, int iStep, int& nBlock,
|
||||
{
|
||||
const RealFunc& pSol = *mySol->getScalarSol();
|
||||
for (j = 1; cit != grid->end_XYZ() && haveXsol; j++, ++cit)
|
||||
field(1,j) = pSol(Vec4(*cit,time));
|
||||
field(1,j) = pSol(Vec4(*cit,time,grid->getParam(j-1)));
|
||||
}
|
||||
else
|
||||
{
|
||||
@ -588,7 +590,7 @@ int SIMoutput::writeGlvS1 (const Vector& psol, int iStep, int& nBlock,
|
||||
cit = grid->begin_XYZ();
|
||||
const RealFunc& sSol = *mySol->getScalarSol();
|
||||
for (j = 1; cit != grid->end_XYZ() && haveXsol; j++, ++cit)
|
||||
field(field.rows(),j) = sSol(Vec4(*cit,time));
|
||||
field(field.rows(),j) = sSol(Vec4(*cit,time,grid->getParam(j-1)));
|
||||
}
|
||||
}
|
||||
|
||||
@ -686,6 +688,9 @@ bool SIMoutput::writeGlvS2 (const Vector& psol, int iStep, int& nBlock,
|
||||
{
|
||||
if (myModel[i]->empty()) continue; // skip empty patches
|
||||
|
||||
if (mySol)
|
||||
mySol->initPatch(myModel[i]->idx);
|
||||
|
||||
if (msgLevel > 1)
|
||||
IFEM::cout <<"Writing secondary solution for patch "<< i+1 << std::endl;
|
||||
|
||||
@ -747,7 +752,7 @@ bool SIMoutput::writeGlvS2 (const Vector& psol, int iStep, int& nBlock,
|
||||
Vec3Vec::const_iterator cit = grid->begin_XYZ();
|
||||
for (j = 1; cit != grid->end_XYZ() && haveAsol; j++, ++cit)
|
||||
{
|
||||
Vec4 Xt(*cit,time);
|
||||
Vec4 Xt(*cit,time,grid->getParam(j-1));
|
||||
if (mySol->hasScalarSol() == 3 || mySol->hasVectorSol() == 3)
|
||||
haveAsol = myProblem->evalSol(lovec,*mySol->getStressSol(),Xt);
|
||||
else if (this->getNoFields() == 1)
|
||||
|
||||
@ -17,6 +17,12 @@
|
||||
#include "Utilities.h"
|
||||
#include "IFEM.h"
|
||||
#include "tinyxml.h"
|
||||
#ifdef HAS_HDF5
|
||||
#include "ProcessAdm.h"
|
||||
#include "HDF5Writer.h"
|
||||
#include "XMLWriter.h"
|
||||
#include "FieldFunctions.h"
|
||||
#endif
|
||||
|
||||
|
||||
AnaSol::AnaSol (RealFunc* s1, VecFunc* s2,
|
||||
@ -107,8 +113,48 @@ template <class T> void parseDerivatives (T* func, const std::string& variables,
|
||||
AnaSol::AnaSol (const TiXmlElement* elem, bool scalarSol)
|
||||
: vecSol(nullptr), vecSecSol(nullptr), stressSol(nullptr)
|
||||
{
|
||||
std::string variables;
|
||||
const char* type = elem->Attribute("type");
|
||||
if (type && !strcasecmp(type,"fields"))
|
||||
parseFieldFunctions(elem, scalarSol);
|
||||
else
|
||||
parseExpressionFunctions(elem, scalarSol);
|
||||
}
|
||||
|
||||
|
||||
AnaSol::~AnaSol ()
|
||||
{
|
||||
for (RealFunc* rf : scalSol)
|
||||
delete rf;
|
||||
for (VecFunc* vf : scalSecSol)
|
||||
delete vf;
|
||||
delete vecSol;
|
||||
delete vecSecSol;
|
||||
delete stressSol;
|
||||
}
|
||||
|
||||
|
||||
void AnaSol::initPatch(size_t pIdx)
|
||||
{
|
||||
for (RealFunc* rf : scalSol)
|
||||
rf->initPatch(pIdx);
|
||||
|
||||
for (VecFunc* rf : scalSecSol)
|
||||
rf->initPatch(pIdx);
|
||||
|
||||
if (vecSol)
|
||||
vecSol->initPatch(pIdx);
|
||||
|
||||
if (vecSecSol)
|
||||
vecSol->initPatch(pIdx);
|
||||
|
||||
if (stressSol)
|
||||
stressSol->initPatch(pIdx);
|
||||
}
|
||||
|
||||
|
||||
void AnaSol::parseExpressionFunctions(const TiXmlElement* elem, bool scalarSol)
|
||||
{
|
||||
std::string variables;
|
||||
const TiXmlElement* var = elem->FirstChildElement("variables");
|
||||
if (var && var->FirstChild())
|
||||
{
|
||||
@ -185,13 +231,110 @@ AnaSol::AnaSol (const TiXmlElement* elem, bool scalarSol)
|
||||
}
|
||||
|
||||
|
||||
AnaSol::~AnaSol ()
|
||||
void AnaSol::parseFieldFunctions(const TiXmlElement* elem, bool scalarSol)
|
||||
{
|
||||
for (RealFunc* rf : scalSol)
|
||||
delete rf;
|
||||
for (VecFunc* vf : scalSecSol)
|
||||
delete vf;
|
||||
delete vecSol;
|
||||
delete vecSecSol;
|
||||
delete stressSol;
|
||||
#ifdef HAS_HDF5
|
||||
std::string file = elem->Attribute("file");
|
||||
int level = 0;
|
||||
utl::getAttribute(elem, "level", level);
|
||||
ProcessAdm adm;
|
||||
XMLWriter xml(file, adm);
|
||||
xml.readInfo();
|
||||
const TiXmlElement* prim = elem->FirstChildElement("primary");
|
||||
if (prim && prim->FirstChild())
|
||||
{
|
||||
std::string primary = prim->FirstChild()->Value();
|
||||
IFEM::cout <<"\tPrimary="<< primary << std::endl;
|
||||
for (const XMLWriter::Entry& it : xml.getEntries()) {
|
||||
if (it.name == primary) {
|
||||
if (scalarSol)
|
||||
scalSol.push_back(new FieldFunction(file, it.basis,
|
||||
primary, it.patches, level));
|
||||
else
|
||||
vecSol = new VecFieldFunction(file, it.basis,
|
||||
primary, it.patches,level);
|
||||
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
prim = elem->FirstChildElement("scalarprimary");
|
||||
if (prim && prim->FirstChild())
|
||||
{
|
||||
std::string primary = prim->FirstChild()->Value();
|
||||
IFEM::cout <<"\tScalar Primary="<< primary << std::endl;
|
||||
for (const XMLWriter::Entry& it : xml.getEntries())
|
||||
if (it.name == primary) {
|
||||
scalSol.push_back(new FieldFunction(file, it.basis,
|
||||
primary, it.patches, level));
|
||||
break;
|
||||
}
|
||||
}
|
||||
const TiXmlElement* sec = elem->FirstChildElement("secondary");
|
||||
if (sec && sec->FirstChild())
|
||||
{
|
||||
std::string secondary = sec->FirstChild()->Value();
|
||||
IFEM::cout <<"\tSecondary="<< secondary << std::endl;
|
||||
// first component runs the show
|
||||
size_t pos = secondary.find_first_of('|');
|
||||
std::string sec;
|
||||
if (pos == std::string::npos)
|
||||
sec = secondary;
|
||||
else
|
||||
sec = secondary.substr(0,pos);
|
||||
for (const XMLWriter::Entry& it : xml.getEntries())
|
||||
if (it.name == sec) {
|
||||
if (scalarSol)
|
||||
scalSecSol.push_back(new VecFieldFunction(file, it.basis, secondary,
|
||||
it.patches, level));
|
||||
else
|
||||
vecSecSol = new TensorFieldFunction(file, it.basis, secondary,
|
||||
it.patches, level);
|
||||
|
||||
break;
|
||||
}
|
||||
}
|
||||
sec = elem->FirstChildElement("scalarsecondary");
|
||||
if (sec && sec->FirstChild())
|
||||
{
|
||||
std::string secondary = sec->FirstChild()->Value();
|
||||
IFEM::cout <<"\tScalar Secondary="<< secondary << std::endl;
|
||||
// first component runs the show
|
||||
size_t pos = secondary.find_first_of('|');
|
||||
std::string sec;
|
||||
if (pos == std::string::npos)
|
||||
sec = secondary;
|
||||
else
|
||||
sec = secondary.substr(0,pos);
|
||||
for (const XMLWriter::Entry& it : xml.getEntries())
|
||||
if (it.name == sec) {
|
||||
scalSecSol.push_back(new VecFieldFunction(file, it.basis, secondary,
|
||||
it.patches, level));
|
||||
break;
|
||||
}
|
||||
}
|
||||
sec = elem->FirstChildElement("stress");
|
||||
if (sec && sec->FirstChild())
|
||||
{
|
||||
std::string secondary = sec->FirstChild()->Value();
|
||||
IFEM::cout <<"\tStress="<< secondary << std::endl;
|
||||
// first component runs the show
|
||||
size_t pos = secondary.find_first_of('|');
|
||||
std::string sec;
|
||||
if (pos == std::string::npos)
|
||||
sec = secondary;
|
||||
else
|
||||
sec = secondary.substr(0,pos);
|
||||
for (const XMLWriter::Entry& it : xml.getEntries())
|
||||
if (it.name == sec) {
|
||||
stressSol = new STensorFieldFunction(file, it.basis, secondary,
|
||||
it.patches, level);
|
||||
|
||||
break;
|
||||
}
|
||||
}
|
||||
#else
|
||||
std::cerr << "AnaSol::parseFieldFunctions(..): Compiled without HDF5 support"
|
||||
<<", no fields read" << std::endl;
|
||||
#endif
|
||||
}
|
||||
|
||||
@ -38,6 +38,12 @@ protected:
|
||||
TensorFunc* vecSecSol; //!< Secondary solution field (vector gradient field)
|
||||
STensorFunc* stressSol; //!< Secondary solution field (stress field)
|
||||
|
||||
//! \brief Parse expression functions from XML definition.
|
||||
void parseExpressionFunctions(const TiXmlElement* elem, bool scalarSol);
|
||||
|
||||
//! \brief Parse field functions from XML definition.
|
||||
void parseFieldFunctions(const TiXmlElement* elem, bool scalarSol);
|
||||
|
||||
public:
|
||||
//! \brief Default constructor initializing all solution fields.
|
||||
//! \param[in] s1 Primary scalar solution field
|
||||
@ -110,6 +116,7 @@ public:
|
||||
//! \brief Returns the scalar solution, if any.
|
||||
RealFunc* getScalarSol(size_t idx = 0) const
|
||||
{ return scalSol.size() <= idx ? nullptr : scalSol[idx]; }
|
||||
|
||||
//! \brief Returns the secondary scalar solution, if any.
|
||||
VecFunc* getScalarSecSol(size_t idx = 0) const
|
||||
{ return scalSecSol.size() <= idx ? nullptr : scalSecSol[idx]; }
|
||||
@ -121,6 +128,9 @@ public:
|
||||
|
||||
//! \brief Returns the stress solution, if any.
|
||||
STensorFunc* getStressSol() const { return stressSol; }
|
||||
|
||||
//! \brief Set patch to use.
|
||||
void initPatch(size_t pIdx);
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
@ -33,6 +33,7 @@ ElementBlock::ElementBlock (size_t nenod)
|
||||
void ElementBlock::resize (size_t nI, size_t nJ, size_t nK)
|
||||
{
|
||||
coord.resize(nI*nJ*nK);
|
||||
params.resize(nI*nJ*nK);
|
||||
if (nen == 2 && nJ < 2 && nK < 2)
|
||||
MMNPC.resize(2*(nI-1));
|
||||
else if (nen == 3 && nK < 2)
|
||||
@ -56,6 +57,7 @@ void ElementBlock::resize (size_t nI, size_t nJ, size_t nK)
|
||||
void ElementBlock::unStructResize (size_t nEl, size_t nPts)
|
||||
{
|
||||
coord.resize(nPts);
|
||||
params.resize(nPts);
|
||||
MMNPC.resize(nen*nEl);
|
||||
MINEX.resize(MMNPC.size()/nen,0);
|
||||
std::iota(MINEX.begin(),MINEX.end(),1);
|
||||
@ -80,6 +82,15 @@ bool ElementBlock::setCoor (size_t i, const Vec3& X)
|
||||
}
|
||||
|
||||
|
||||
bool ElementBlock::setParams (size_t i, Real u, Real v, Real w)
|
||||
{
|
||||
if (i >= params.size()) return false;
|
||||
|
||||
params[i] = {u,v,w};
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
bool ElementBlock::setNode (size_t i, int nodeNumb)
|
||||
{
|
||||
if (i >= MMNPC.size()) return false;
|
||||
|
||||
@ -15,8 +15,10 @@
|
||||
#define _ELEMENTBLOCK_H
|
||||
|
||||
#include "Vec3.h"
|
||||
#include <array>
|
||||
#include <vector>
|
||||
|
||||
|
||||
/*!
|
||||
\brief Class for storage of a standard FE grid block.
|
||||
*/
|
||||
@ -45,6 +47,8 @@ public:
|
||||
//! \brief Defines the coordinates of node \a i.
|
||||
bool setCoor(size_t i, Real x, Real y, Real z)
|
||||
{ return this->setCoor(i,Vec3(x,y,z)); }
|
||||
//! \brief Defines the parameter values for node \a i.
|
||||
bool setParams(size_t i, Real u, Real v, Real w = Real(0));
|
||||
|
||||
//! \brief Defines the global number of element node \a i.
|
||||
bool setNode(size_t i, int nodeNumb);
|
||||
@ -76,14 +80,18 @@ public:
|
||||
//! \brief Returns the coordinate of a given node.
|
||||
const Vec3& getCoord(size_t i) const { return coord[i]; }
|
||||
|
||||
//! \brief Returns a pointer to the parameter values of a given node.
|
||||
const Real* getParam(size_t i) const { return params[i].data(); }
|
||||
|
||||
//! \brief Returns a pointer to the element connectivity array.
|
||||
const int* getElements() const { return &MMNPC.front(); }
|
||||
const int* getElements() const { return MMNPC.data(); }
|
||||
|
||||
private:
|
||||
std::vector<Vec3> coord; //!< Vector of nodal coordinates
|
||||
std::vector<int> MMNPC; //!< Matrix of Matrices of Nodal Point Correspondance
|
||||
std::vector<int> MINEX; //!< Matrix of Internal to External element numbers
|
||||
size_t nen; //!< Number of Element Nodes
|
||||
std::vector<std::array<Real,3>> params; //!< Parameter values where relevant
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
@ -14,9 +14,12 @@
|
||||
#include "FieldFunctions.h"
|
||||
#include "ASMbase.h"
|
||||
#include "ASM2D.h"
|
||||
#include "ASM3D.h"
|
||||
#include "Field.h"
|
||||
#include "Fields.h"
|
||||
#include "ProcessAdm.h"
|
||||
#include "Vec3.h"
|
||||
#include "StringUtils.h"
|
||||
#ifdef HAS_HDF5
|
||||
#include "HDF5Writer.h"
|
||||
#endif
|
||||
@ -25,46 +28,197 @@
|
||||
|
||||
FieldFunction::FieldFunction (const std::string& fileName,
|
||||
const std::string& basisName,
|
||||
const std::string& fieldName)
|
||||
const std::string& fieldName,
|
||||
size_t npatches, int level)
|
||||
{
|
||||
pidx = 0;
|
||||
#ifdef HAS_HDF5
|
||||
HDF5Writer hdf5(fileName,ProcessAdm(),true,true);
|
||||
ProcessAdm adm;
|
||||
HDF5Writer hdf5(fileName,adm,true,true);
|
||||
|
||||
std::string g2;
|
||||
std::stringstream str;
|
||||
hdf5.readString("0/basis/"+basisName+"/1",g2);
|
||||
str << g2;
|
||||
for (size_t p = 0; p < npatches; ++p) {
|
||||
std::string g2;
|
||||
std::stringstream str;
|
||||
std::stringstream str2;
|
||||
str2 << level << "/basis/" << basisName << "/" << p+1;
|
||||
hdf5.readString(str2.str(),g2);
|
||||
str << g2;
|
||||
std::string head = g2.substr(0,9);
|
||||
std::string head2 = g2.substr(0,12);
|
||||
ASMbase* pch = nullptr;
|
||||
if (head == "200 1 0 0")
|
||||
pch = ASM2D::create(ASM::Spline);
|
||||
else if (head == "700 1 0 0")
|
||||
pch = ASM3D::create(ASM::Spline);
|
||||
else if (head2 == "# LRSPLINE S")
|
||||
pch = ASM2D::create(ASM::LRSpline);
|
||||
else if (head2 == "# LRSPLINE V")
|
||||
pch = ASM3D::create(ASM::LRSpline);
|
||||
|
||||
pch = ASM2D::create(ASM::Spline);
|
||||
pch->read(str);
|
||||
|
||||
Vector coefs;
|
||||
hdf5.readVector(0,fieldName,1,coefs);
|
||||
field = Field::create(pch,coefs);
|
||||
if (pch)
|
||||
{
|
||||
pch->read(str);
|
||||
Vector coefs;
|
||||
hdf5.readVector(level,fieldName,p+1,coefs);
|
||||
field.push_back(Field::create(pch,coefs));
|
||||
patch.push_back(pch);
|
||||
}
|
||||
else {
|
||||
std::cerr <<" *** FieldFunction: Unknown basis type, no function created."
|
||||
<< std::endl;
|
||||
return;
|
||||
}
|
||||
}
|
||||
#else
|
||||
std::cerr <<"WARNING: Compiled without HDF5 support,"
|
||||
<<" field function is not instanciated."<< std::endl;
|
||||
field = nullptr;
|
||||
pch = nullptr;
|
||||
<<" field function is not instantiated."<< std::endl;
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
FieldFunction::~FieldFunction ()
|
||||
{
|
||||
delete field;
|
||||
delete pch;
|
||||
for (Field* f : field)
|
||||
delete f;
|
||||
for (ASMbase* pch : patch)
|
||||
delete pch;
|
||||
}
|
||||
|
||||
|
||||
Real FieldFunction::evaluate (const Vec3& X) const
|
||||
{
|
||||
if (!field)
|
||||
if (pidx >= field.size() || !field[pidx])
|
||||
return Real(0);
|
||||
|
||||
const Vec4* x4 = dynamic_cast<const Vec4*>(&X);
|
||||
if (x4 && x4->idx > 0)
|
||||
return field->valueNode(x4->idx);
|
||||
if (x4)
|
||||
return x4->idx > 0 ? field[pidx]->valueNode(x4->idx)
|
||||
: field[pidx]->valueCoor(*x4);
|
||||
|
||||
return field->valueCoor(X);
|
||||
return field[pidx]->valueCoor(X);
|
||||
}
|
||||
|
||||
|
||||
VecFieldFuncBase::VecFieldFuncBase (const std::string& fileName,
|
||||
const std::string& basisName,
|
||||
const std::string& fieldName,
|
||||
size_t npatches, int level)
|
||||
{
|
||||
pidx = 0;
|
||||
#ifdef HAS_HDF5
|
||||
HDF5Writer hdf5(fileName,ProcessAdm(),true,true);
|
||||
|
||||
std::vector<std::string> fieldNames = splitString(fieldName,
|
||||
[](int c) { return c == '|' ? 1 : 0; });
|
||||
|
||||
for (size_t p = 0; p < npatches; ++p) {
|
||||
std::string g2;
|
||||
std::stringstream str;
|
||||
std::stringstream str2;
|
||||
str2 << level << "/basis/" << basisName << "/" << p+1;
|
||||
hdf5.readString(str2.str(),g2);
|
||||
str << g2;
|
||||
std::string head = g2.substr(0,9);
|
||||
std::string head2 = g2.substr(0,12);
|
||||
ASMbase* pch = nullptr;
|
||||
if (head == "200 1 0 0")
|
||||
pch = ASM2D::create(ASM::Spline, std::max(fieldNames.size(), 2LU));
|
||||
else if (head == "700 1 0 0")
|
||||
pch = ASM3D::create(ASM::Spline, std::max(fieldNames.size(), 3LU));
|
||||
else if (head2 == "# LRSPLINE S")
|
||||
pch = ASM2D::create(ASM::LRSpline, std::max(fieldNames.size(), 2LU));
|
||||
else if (head2 == "# LRSPLINE V")
|
||||
pch = ASM3D::create(ASM::LRSpline, std::max(fieldNames.size(), 3LU));
|
||||
|
||||
if (pch)
|
||||
{
|
||||
pch->read(str);
|
||||
Vector coefs;
|
||||
if (fieldNames.size() == 1)
|
||||
hdf5.readVector(level,fieldName,p+1,coefs);
|
||||
else {
|
||||
Vectors coef(fieldNames.size());
|
||||
for (size_t i = 0; i < fieldNames.size(); ++i)
|
||||
hdf5.readVector(level,fieldNames[i],p+1,coef[i]);
|
||||
coefs.reserve(coef.size()*coef.front().size());
|
||||
for (size_t i = 0; i < coef.front().size(); ++i)
|
||||
for (size_t j = 0; j < coef.size(); ++j)
|
||||
coefs.push_back(coef[j][i]);
|
||||
}
|
||||
field.push_back(Fields::create(pch,coefs,1,pch->getNoFields(1)));
|
||||
patch.push_back(pch);
|
||||
}
|
||||
else {
|
||||
std::cerr <<" *** VecFieldFunction: Unknown basis type, no function created."
|
||||
<< std::endl;
|
||||
return;
|
||||
}
|
||||
}
|
||||
#else
|
||||
std::cerr <<"WARNING: Compiled without HDF5 support,"
|
||||
<<" vector field function is not instantiated."<< std::endl;
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
VecFieldFuncBase::~VecFieldFuncBase ()
|
||||
{
|
||||
for (Fields* f : field)
|
||||
delete f;
|
||||
for (ASMbase* pch : patch)
|
||||
delete pch;
|
||||
}
|
||||
|
||||
|
||||
Vec3 VecFieldFunction::evaluate (const Vec3& X) const
|
||||
{
|
||||
if (pidx >= field.size() || !field[pidx])
|
||||
return Vec3();
|
||||
|
||||
Vector vals;
|
||||
const Vec4* x4 = dynamic_cast<const Vec4*>(&X);
|
||||
if (!x4)
|
||||
field[pidx]->valueCoor(X, vals);
|
||||
else if (x4->idx > 0)
|
||||
field[pidx]->valueNode(x4->idx, vals);
|
||||
else
|
||||
field[pidx]->valueCoor(*x4, vals);
|
||||
|
||||
return Vec3(vals.ptr(), patch[pidx]->getNoSpaceDim());
|
||||
}
|
||||
|
||||
|
||||
Tensor TensorFieldFunction::evaluate (const Vec3& X) const
|
||||
{
|
||||
if (pidx >= field.size() || !field[pidx])
|
||||
return Tensor(3);
|
||||
|
||||
Vector vals;
|
||||
const Vec4* x4 = dynamic_cast<const Vec4*>(&X);
|
||||
if (!x4)
|
||||
field[pidx]->valueCoor(X, vals);
|
||||
else if (x4->idx > 0)
|
||||
field[pidx]->valueNode(x4->idx, vals);
|
||||
else
|
||||
field[pidx]->valueCoor(*x4, vals);
|
||||
|
||||
return vals;
|
||||
}
|
||||
|
||||
|
||||
SymmTensor STensorFieldFunction::evaluate (const Vec3& X) const
|
||||
{
|
||||
if (pidx >= field.size() || !field[pidx])
|
||||
return SymmTensor(3);
|
||||
|
||||
Vector vals;
|
||||
const Vec4* x4 = dynamic_cast<const Vec4*>(&X);
|
||||
if (!x4)
|
||||
field[pidx]->valueCoor(X, vals);
|
||||
else if (x4->idx > 0)
|
||||
field[pidx]->valueNode(x4->idx, vals);
|
||||
else
|
||||
field[pidx]->valueCoor(*x4, vals);
|
||||
|
||||
return vals;
|
||||
}
|
||||
|
||||
@ -15,9 +15,11 @@
|
||||
#define _FIELD_FUNCTIONS_H
|
||||
|
||||
#include "Function.h"
|
||||
#include "TensorFunction.h"
|
||||
#include <string>
|
||||
|
||||
class Field;
|
||||
class Fields;
|
||||
class ASMbase;
|
||||
|
||||
|
||||
@ -27,25 +29,165 @@ class ASMbase;
|
||||
|
||||
class FieldFunction : public RealFunc
|
||||
{
|
||||
Field* field; //!< Pointer to the scalar field to be evaluated.
|
||||
ASMbase* pch; //!< Pointer to the patch on which the field is defined over
|
||||
std::vector<Field*> field; //!< Pointer to the scalar field to be evaluated
|
||||
std::vector<ASMbase*> patch; //!< Pointer to the patch on which the field is defined
|
||||
|
||||
public:
|
||||
//! \brief Default constructor.
|
||||
explicit FieldFunction(Field* f = nullptr) : field(f), pch(nullptr) {}
|
||||
explicit FieldFunction(Field* f = nullptr) : field{f}, pidx(0) {}
|
||||
//! \brief Constructor creating a field from a provided HDF5 file.
|
||||
//! \param[in] fileName Name of the HDF5-file
|
||||
//! \param[in] basisName Name of the basis which the field values refer to
|
||||
//! \param[in] fieldName Name of the field in the HDF5-file
|
||||
//! \param[in] npatches Number of patches to read for
|
||||
//! \param[in] level Level to read for
|
||||
FieldFunction(const std::string& fileName,
|
||||
const std::string& basisName,
|
||||
const std::string& fieldName);
|
||||
//! \brief The destructor deletes the scalar field.
|
||||
const std::string& fieldName,
|
||||
size_t npatches = 1, int level = 0);
|
||||
|
||||
//! \brief The destructor deletes the scalar fields.
|
||||
virtual ~FieldFunction();
|
||||
|
||||
//! \brief Set currently active patch.
|
||||
virtual void initPatch(size_t pIdx) { if (pIdx < field.size()) pidx = pIdx; }
|
||||
|
||||
protected:
|
||||
//! \brief Evaluates the scalar field function.
|
||||
virtual Real evaluate(const Vec3& X) const;
|
||||
|
||||
size_t pidx; //!< Current patch index
|
||||
};
|
||||
|
||||
/*!
|
||||
\brief Base class for multi-value spatial function, defined trough a vector field.
|
||||
*/
|
||||
|
||||
class VecFieldFuncBase
|
||||
{
|
||||
protected:
|
||||
//! \brief Default constructor.
|
||||
VecFieldFuncBase(Fields* f = nullptr) : field{f} {}
|
||||
|
||||
//! \brief Constructor creating a field from a provided HDF5 file.
|
||||
//! \param[in] fileName Name of the HDF5-file
|
||||
//! \param[in] basisName Name of the basis which the field values refer to
|
||||
//! \param[in] fieldName Name of the field in the HDF5-file
|
||||
//! \param[in] patch Number of patches to read for
|
||||
//! \param[in] level Level to read for
|
||||
VecFieldFuncBase(const std::string& fileName,
|
||||
const std::string& basisName,
|
||||
const std::string& fieldName,
|
||||
size_t npatches = 1, int level = 0);
|
||||
|
||||
//! \brief The destructor deletes the vector field.
|
||||
virtual ~VecFieldFuncBase();
|
||||
|
||||
std::vector<Fields*> field; //!< Pointer to the vector field to be evaluated
|
||||
std::vector<ASMbase*> patch; //!< Pointer to the patch on which the field is defined
|
||||
size_t pidx; //!< Current patch index
|
||||
};
|
||||
|
||||
|
||||
/*!
|
||||
\brief A vector-valued spatial function, defined trough a vector field.
|
||||
*/
|
||||
|
||||
class VecFieldFunction : public VecFunc, public VecFieldFuncBase
|
||||
{
|
||||
public:
|
||||
//! \brief Default constructor.
|
||||
explicit VecFieldFunction(Fields* f = nullptr) : VecFieldFuncBase(f) {}
|
||||
|
||||
//! \brief Constructor creating a field from a provided HDF5 file.
|
||||
//! \param[in] fileName Name of the HDF5-file
|
||||
//! \param[in] basisName Name of the basis which the field values refer to
|
||||
//! \param[in] fieldName Name of the field in the HDF5-file
|
||||
//! \param[in] nPatches Number of patches to read for
|
||||
//! \param[in] level Level to read for
|
||||
VecFieldFunction(const std::string& fileName,
|
||||
const std::string& basisName,
|
||||
const std::string& fieldName,
|
||||
size_t nPatches = 1, int level = 0) :
|
||||
VecFieldFuncBase(fileName, basisName, fieldName, nPatches, level) {}
|
||||
|
||||
//! \brief Empty destructor.
|
||||
virtual ~VecFieldFunction() {}
|
||||
|
||||
//! \brief Set currently active patch.
|
||||
virtual void initPatch(size_t pIdx) { if (pIdx < field.size()) pidx = pIdx; }
|
||||
|
||||
protected:
|
||||
//! \brief Evaluates the vectorial field function.
|
||||
virtual Vec3 evaluate(const Vec3& X) const;
|
||||
};
|
||||
|
||||
|
||||
/*!
|
||||
\brief A tensor-valued spatial function, defined trough a vector field.
|
||||
*/
|
||||
|
||||
class TensorFieldFunction : public TensorFunc, public VecFieldFuncBase
|
||||
{
|
||||
public:
|
||||
//! \brief Default constructor.
|
||||
explicit TensorFieldFunction(Fields* f = nullptr) : VecFieldFuncBase(f) {}
|
||||
|
||||
//! \brief Constructor creating a field from a provided HDF5 file.
|
||||
//! \param[in] fileName Name of the HDF5-file
|
||||
//! \param[in] basisName Name of the basis which the field values refer to
|
||||
//! \param[in] fieldName Name of the field in the HDF5-file
|
||||
//! \param[in] nPatches Number of patches to read for
|
||||
//! \param[in] level Level to read for
|
||||
TensorFieldFunction(const std::string& fileName,
|
||||
const std::string& basisName,
|
||||
const std::string& fieldName,
|
||||
size_t nPatches = 1, int level = 0) :
|
||||
VecFieldFuncBase(fileName, basisName, fieldName, nPatches, level) {}
|
||||
|
||||
//! \brief Empty destructor.
|
||||
virtual ~TensorFieldFunction() {}
|
||||
|
||||
//! \brief Set currently active patch.
|
||||
virtual void initPatch(size_t pIdx) { if (pIdx < field.size()) pidx = pIdx; }
|
||||
|
||||
protected:
|
||||
//! \brief Evaluates the tensorial field function.
|
||||
virtual Tensor evaluate(const Vec3& X) const;
|
||||
};
|
||||
|
||||
|
||||
/*!
|
||||
\brief A symmtensor-valued spatial function, defined trough a vector field.
|
||||
*/
|
||||
|
||||
class STensorFieldFunction : public STensorFunc, public VecFieldFuncBase
|
||||
{
|
||||
public:
|
||||
//! \brief Default constructor.
|
||||
explicit STensorFieldFunction(Fields* f = nullptr) : VecFieldFuncBase(f) {}
|
||||
|
||||
//! \brief Constructor creating a field from a provided HDF5 file.
|
||||
//! \param[in] fileName Name of the HDF5-file
|
||||
//! \param[in] basisName Name of the basis which the field values refer to
|
||||
//! \param[in] fieldName Name of the field in the HDF5-file
|
||||
//! \param[in] nPatches Number of patches to read for
|
||||
//! \param[in] level Level to read for
|
||||
STensorFieldFunction(const std::string& fileName,
|
||||
const std::string& basisName,
|
||||
const std::string& fieldName,
|
||||
size_t nPatches = 1, int level = 0) :
|
||||
VecFieldFuncBase(fileName, basisName, fieldName, nPatches, level) {}
|
||||
|
||||
//! \brief Empty destructor.
|
||||
virtual ~STensorFieldFunction() {}
|
||||
|
||||
//! \brief Set currently active patch.
|
||||
virtual void initPatch(size_t pIdx) { if (pIdx < field.size()) pidx = pIdx; }
|
||||
|
||||
protected:
|
||||
//! \brief Evaluates the tensorial field function.
|
||||
virtual SymmTensor evaluate(const Vec3& X) const;
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
@ -142,6 +142,9 @@ public:
|
||||
//! \brief Returns the number of components of the return value.
|
||||
size_t dim() const { return ncmp; }
|
||||
|
||||
//! \brief Set currently active patch.
|
||||
virtual void initPatch(size_t pIdx) {}
|
||||
|
||||
protected:
|
||||
size_t ncmp; //!< Number of components in the return value
|
||||
};
|
||||
|
||||
187
src/Utility/Test/TestFieldFunctions.C
Normal file
187
src/Utility/Test/TestFieldFunctions.C
Normal file
@ -0,0 +1,187 @@
|
||||
//==============================================================================
|
||||
//!
|
||||
//! \file TestFieldFunctions.C
|
||||
//!
|
||||
//! \date Nov 7 2017
|
||||
//!
|
||||
//! \author Arne Morten Kvarving / SINTEF
|
||||
//!
|
||||
//! \brief Tests for parsing of field functions.
|
||||
//!
|
||||
//==============================================================================
|
||||
|
||||
#include "FieldFunctions.h"
|
||||
|
||||
#include "gtest/gtest.h"
|
||||
|
||||
|
||||
TEST(TestFieldFunctions, 2D1P)
|
||||
{
|
||||
FieldFunction f2D_scalar("src/Utility/Test/refdata/Field2D-1P",
|
||||
"Stokes-2", "v", 1, 0);
|
||||
|
||||
VecFieldFunction f2D_vec("src/Utility/Test/refdata/Field2D-1P",
|
||||
"Stokes-1", "u", 1, 0);
|
||||
|
||||
TensorFieldFunction f2D_ten("src/Utility/Test/refdata/Field2D-1P",
|
||||
"Stokes-1", "u_x,x|u_x,y|u_y,x|u_y,y", 1, 0);
|
||||
|
||||
STensorFieldFunction f2D_sten("src/Utility/Test/refdata/Field2D-1P",
|
||||
"Stokes-1", "u_x,x|u_y,y|u_x,y", 1, 0);
|
||||
|
||||
double param[3] = {0.5, 0.5, 0.0};
|
||||
Vec4 X(param);
|
||||
double scal = f2D_scalar(X);
|
||||
double x = 1.0, y = 1.0;
|
||||
EXPECT_NEAR(scal, x, 1e-14);
|
||||
Vec3 vec = f2D_vec(X);
|
||||
EXPECT_NEAR(vec[0], x, 1e-14);
|
||||
EXPECT_NEAR(vec[1], -y, 1e-14);
|
||||
Tensor ten = f2D_ten(X); // see above
|
||||
EXPECT_NEAR(ten(1,1), 1.0, 1e-14);
|
||||
EXPECT_NEAR(ten(1,2), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(2,1), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(2,2), -1.0, 1e-14);
|
||||
SymmTensor sten = f2D_sten(X);
|
||||
EXPECT_NEAR(sten(1,1), 1.0, 1e-14);
|
||||
EXPECT_NEAR(sten(1,2), 0.0, 1e-14);
|
||||
EXPECT_NEAR(sten(2,1), 0.0, 1e-14);
|
||||
EXPECT_NEAR(sten(2,2), -1.0, 1e-14);
|
||||
|
||||
param[0] = 1.0;
|
||||
param[1] = 0.25;
|
||||
scal = f2D_scalar(X);
|
||||
x = 2.0, y = 0.5;
|
||||
vec = f2D_vec(X);
|
||||
EXPECT_NEAR(scal, x, 1e-14);
|
||||
EXPECT_NEAR(vec[0], x, 1e-14);
|
||||
EXPECT_NEAR(vec[1], -y, 1e-14);
|
||||
}
|
||||
|
||||
|
||||
TEST(TestFieldFunctions, 2D2P)
|
||||
{
|
||||
FieldFunction f2D_scalar("src/Utility/Test/refdata/Field2D-2P",
|
||||
"Stokes-2", "v", 2, 0);
|
||||
|
||||
VecFieldFunction f2D_vec("src/Utility/Test/refdata/Field2D-2P",
|
||||
"Stokes-1", "u", 2, 0);
|
||||
|
||||
double param[3] = {0.5, 0.5, 0.0};
|
||||
Vec4 X(param);
|
||||
for (size_t pIdx = 0; pIdx < 2; ++pIdx) {
|
||||
f2D_scalar.initPatch(pIdx);
|
||||
double scal = f2D_scalar(X);
|
||||
EXPECT_NEAR(scal, 0.5 + pIdx, 1e-14);
|
||||
f2D_vec.initPatch(pIdx);
|
||||
Vec3 vec = f2D_vec(X);
|
||||
EXPECT_NEAR(vec[0], 0.5 + pIdx, 1e-14);
|
||||
EXPECT_NEAR(vec[1], -1.0, 1e-14);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
TEST(TestFieldFunctions, 3D1P)
|
||||
{
|
||||
FieldFunction f3D_scalar("src/Utility/Test/refdata/Field3D-1P",
|
||||
"Stokes-2", "v", 1, 0);
|
||||
|
||||
VecFieldFunction f3D_vec("src/Utility/Test/refdata/Field3D-1P",
|
||||
"Stokes-1", "u", 1, 0);
|
||||
|
||||
TensorFieldFunction f3D_ten("src/Utility/Test/refdata/Field3D-1P",
|
||||
"Stokes-1", "u_x,x|u_y,x|u_z,x|u_x,y|u_y,y|u_z,y|u_x,z|u_x,y|u_z,z", 1, 0);
|
||||
|
||||
STensorFieldFunction f3D_sten("src/Utility/Test/refdata/Field3D-1P",
|
||||
"Stokes-1", "u_x,x|u_y,y|u_z,z|u_x,y|u_z,x|u_x,y", 1, 0);
|
||||
|
||||
double param[3] = {0.25, 0.25, 0.25};
|
||||
Vec4 X(param);
|
||||
double scal = f3D_scalar(X);
|
||||
double x = 0.5;
|
||||
double y = 0.5;
|
||||
double z = 0.5;
|
||||
EXPECT_NEAR(scal, 0.0, 1e-14);
|
||||
Vec3 vec = f3D_vec(X);
|
||||
EXPECT_NEAR(vec[0], y*(1.0-y), 1e-14); // y*(1-y)
|
||||
EXPECT_NEAR(vec[1], z*(1.0-z), 1e-14); // z*(1-z)
|
||||
EXPECT_NEAR(vec[2], x*(1.0-x), 1e-14); // x*(1-x)
|
||||
param[0] = param[1] = param[2] = 0.5;
|
||||
Tensor ten = f3D_ten(X);
|
||||
x = y = z = 1.0;
|
||||
EXPECT_NEAR(ten(1,1), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(1,2), 1.0-2*y, 1e-14);
|
||||
EXPECT_NEAR(ten(1,3), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(2,1), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(2,2), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(2,3), 1.0-2*z, 1e-14);
|
||||
EXPECT_NEAR(ten(3,1), 1.0-2*x, 1e-14);
|
||||
EXPECT_NEAR(ten(3,2), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(3,3), 0.0, 1e-14);
|
||||
|
||||
SymmTensor sten = f3D_sten(X);
|
||||
EXPECT_NEAR(sten(1,1), 0.0, 1e-14);
|
||||
EXPECT_NEAR(sten(1,2), 1.0-2*x, 1e-14);
|
||||
EXPECT_NEAR(sten(1,3), 1.0-2*y, 1e-14);
|
||||
EXPECT_NEAR(sten(2,1), 1.0-2*x, 1e-14);
|
||||
EXPECT_NEAR(sten(2,2), 0.0, 1e-14);
|
||||
EXPECT_NEAR(sten(2,3), 1.0-2*x, 1e-14);
|
||||
EXPECT_NEAR(sten(3,1), 1.0-2*y, 1e-14);
|
||||
EXPECT_NEAR(sten(3,2), 1.0-2*x, 1e-14);
|
||||
EXPECT_NEAR(sten(3,3), 0.0, 1e-14);
|
||||
}
|
||||
|
||||
|
||||
TEST(TestFieldFunctions, 3D2P)
|
||||
{
|
||||
FieldFunction f3D_scalar("src/Utility/Test/refdata/Field3D-2P",
|
||||
"Stokes-2", "v", 2, 0);
|
||||
|
||||
VecFieldFunction f3D_vec("src/Utility/Test/refdata/Field3D-2P",
|
||||
"Stokes-1", "u", 2, 0);
|
||||
|
||||
TensorFieldFunction f3D_ten("src/Utility/Test/refdata/Field3D-2P",
|
||||
"Stokes-1", "u_x,x|u_y,x|u_z,x|u_x,y|u_y,y|u_z,y|u_x,z|u_x,y|u_z,z", 2, 0);
|
||||
|
||||
STensorFieldFunction f3D_sten("src/Utility/Test/refdata/Field3D-2P",
|
||||
"Stokes-1", "u_x,x|u_y,y|u_z,z|u_x,y|u_y,z|u_z,x", 2, 0);
|
||||
|
||||
double param[3] = {0.25, 0.25, 0.25};
|
||||
Vec4 X(param);
|
||||
for (size_t pIdx = 0; pIdx < 2; ++pIdx) {
|
||||
f3D_scalar.initPatch(pIdx);
|
||||
double scal = f3D_scalar(X);
|
||||
double x = 0.25 + pIdx;
|
||||
double y = 0.5;
|
||||
double z = 0.5;
|
||||
EXPECT_NEAR(scal, 1.0-2*x, 1e-14);
|
||||
f3D_vec.initPatch(pIdx);
|
||||
Vec3 vec = f3D_vec(X);
|
||||
EXPECT_NEAR(vec[0], y*(1.0-y), 1e-14); // y*(1-y)
|
||||
EXPECT_NEAR(vec[1], z*(1.0-z), 1e-14); // z*(1-z)
|
||||
EXPECT_NEAR(vec[2], x*(1.0-x), 1e-14); // x*(1-x)
|
||||
f3D_ten.initPatch(pIdx);
|
||||
Tensor ten = f3D_ten(X);
|
||||
EXPECT_NEAR(ten(1,1), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(1,2), 1.0-2*y, 1e-14);
|
||||
EXPECT_NEAR(ten(1,3), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(2,1), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(2,2), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(2,3), 1.0-2*z, 1e-14);
|
||||
EXPECT_NEAR(ten(3,1), 1.0-2*x, 1e-14);
|
||||
EXPECT_NEAR(ten(3,2), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(3,3), 0.0, 1e-14);
|
||||
|
||||
f3D_sten.initPatch(pIdx);
|
||||
SymmTensor sten = f3D_sten(X);
|
||||
EXPECT_NEAR(sten(1,1), 0.0, 1e-14);
|
||||
EXPECT_NEAR(sten(1,2), 1.0-2*y, 1e-14);
|
||||
EXPECT_NEAR(sten(1,3), 1.0-2*x, 1e-14);
|
||||
EXPECT_NEAR(sten(2,1), 1.0-2*y, 1e-14);
|
||||
EXPECT_NEAR(sten(2,2), 0.0, 1e-14);
|
||||
EXPECT_NEAR(sten(2,3), 1.0-2*z, 1e-14);
|
||||
EXPECT_NEAR(sten(3,1), 1.0-2*x, 1e-14);
|
||||
EXPECT_NEAR(sten(3,2), 1.0-2*z, 1e-14);
|
||||
EXPECT_NEAR(sten(3,3), 0.0, 1e-14);
|
||||
}
|
||||
}
|
||||
187
src/Utility/Test/TestFieldFunctionsLR.C
Normal file
187
src/Utility/Test/TestFieldFunctionsLR.C
Normal file
@ -0,0 +1,187 @@
|
||||
//==============================================================================
|
||||
//!
|
||||
//! \file TestFieldFunctionsLR.C
|
||||
//!
|
||||
//! \date Nov 7 2017
|
||||
//!
|
||||
//! \author Arne Morten Kvarving / SINTEF
|
||||
//!
|
||||
//! \brief Tests for parsing of LR field functions.
|
||||
//!
|
||||
//==============================================================================
|
||||
|
||||
#include "FieldFunctions.h"
|
||||
|
||||
#include "gtest/gtest.h"
|
||||
|
||||
|
||||
TEST(TestFieldFunctions, 2D1PLR)
|
||||
{
|
||||
FieldFunction f2D_scalar("src/Utility/Test/refdata/Field2D-1PLR",
|
||||
"Stokes-2", "v", 1, 0);
|
||||
|
||||
VecFieldFunction f2D_vec("src/Utility/Test/refdata/Field2D-1PLR",
|
||||
"Stokes-1", "u", 1, 0);
|
||||
|
||||
TensorFieldFunction f2D_ten("src/Utility/Test/refdata/Field2D-1PLR",
|
||||
"Stokes-1", "u_x,x|u_x,y|u_y,x|u_y,y", 1, 0);
|
||||
|
||||
STensorFieldFunction f2D_sten("src/Utility/Test/refdata/Field2D-1PLR",
|
||||
"Stokes-1", "u_x,x|u_y,y|u_x,y", 1, 0);
|
||||
|
||||
double param[3] = {0.5, 0.5, 0.0};
|
||||
Vec4 X(param);
|
||||
double scal = f2D_scalar(X);
|
||||
double x = 1.0, y = 1.0;
|
||||
EXPECT_NEAR(scal, x, 1e-14);
|
||||
Vec3 vec = f2D_vec(X);
|
||||
EXPECT_NEAR(vec[0], x, 1e-14);
|
||||
EXPECT_NEAR(vec[1], -y, 1e-14);
|
||||
Tensor ten = f2D_ten(X); // see above
|
||||
EXPECT_NEAR(ten(1,1), 1.0, 1e-14);
|
||||
EXPECT_NEAR(ten(1,2), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(2,1), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(2,2), -1.0, 1e-14);
|
||||
SymmTensor sten = f2D_sten(X);
|
||||
EXPECT_NEAR(sten(1,1), 1.0, 1e-14);
|
||||
EXPECT_NEAR(sten(1,2), 0.0, 1e-14);
|
||||
EXPECT_NEAR(sten(2,1), 0.0, 1e-14);
|
||||
EXPECT_NEAR(sten(2,2), -1.0, 1e-14);
|
||||
|
||||
param[0] = 1.0;
|
||||
param[1] = 0.25;
|
||||
scal = f2D_scalar(X);
|
||||
x = 2.0, y = 0.5;
|
||||
vec = f2D_vec(X);
|
||||
EXPECT_NEAR(scal, x, 1e-14);
|
||||
EXPECT_NEAR(vec[0], x, 1e-14);
|
||||
EXPECT_NEAR(vec[1], -y, 1e-14);
|
||||
}
|
||||
|
||||
|
||||
TEST(TestFieldFunctions, 2D2PLR)
|
||||
{
|
||||
FieldFunction f2D_scalar("src/Utility/Test/refdata/Field2D-2PLR",
|
||||
"Stokes-2", "v", 2, 0);
|
||||
|
||||
VecFieldFunction f2D_vec("src/Utility/Test/refdata/Field2D-2PLR",
|
||||
"Stokes-1", "u", 2, 0);
|
||||
|
||||
double param[3] = {0.5, 0.5, 0.0};
|
||||
Vec4 X(param);
|
||||
for (size_t pIdx = 0; pIdx < 2; ++pIdx) {
|
||||
f2D_scalar.initPatch(pIdx);
|
||||
double scal = f2D_scalar(X);
|
||||
EXPECT_NEAR(scal, 0.5 + pIdx, 1e-14);
|
||||
f2D_vec.initPatch(pIdx);
|
||||
Vec3 vec = f2D_vec(X);
|
||||
EXPECT_NEAR(vec[0], 0.5 + pIdx, 1e-14);
|
||||
EXPECT_NEAR(vec[1], -1.0, 1e-14);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
TEST(TestFieldFunctions, 3D1PLR)
|
||||
{
|
||||
FieldFunction f3D_scalar("src/Utility/Test/refdata/Field3D-1PLR",
|
||||
"Stokes-2", "v", 1, 0);
|
||||
|
||||
VecFieldFunction f3D_vec("src/Utility/Test/refdata/Field3D-1PLR",
|
||||
"Stokes-1", "u", 1, 0);
|
||||
|
||||
TensorFieldFunction f3D_ten("src/Utility/Test/refdata/Field3D-1PLR",
|
||||
"Stokes-1", "u_x,x|u_y,x|u_z,x|u_x,y|u_y,y|u_z,y|u_x,z|u_x,y|u_z,z", 1, 0);
|
||||
|
||||
STensorFieldFunction f3D_sten("src/Utility/Test/refdata/Field3D-1PLR",
|
||||
"Stokes-1", "u_x,x|u_y,y|u_z,z|u_x,y|u_z,x|u_x,y", 1, 0);
|
||||
|
||||
double param[3] = {0.25, 0.25, 0.25};
|
||||
Vec4 X(param);
|
||||
double scal = f3D_scalar(X);
|
||||
double x = 0.5;
|
||||
double y = 0.5;
|
||||
double z = 0.5;
|
||||
EXPECT_NEAR(scal, 0.0, 1e-14);
|
||||
Vec3 vec = f3D_vec(X);
|
||||
EXPECT_NEAR(vec[0], y*(1.0-y), 1e-14); // y*(1-y)
|
||||
EXPECT_NEAR(vec[1], z*(1.0-z), 1e-14); // z*(1-z)
|
||||
EXPECT_NEAR(vec[2], x*(1.0-x), 1e-14); // x*(1-x)
|
||||
param[0] = param[1] = param[2] = 0.5;
|
||||
Tensor ten = f3D_ten(X);
|
||||
x = y = z = 1.0;
|
||||
EXPECT_NEAR(ten(1,1), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(1,2), 1.0-2*y, 1e-14);
|
||||
EXPECT_NEAR(ten(1,3), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(2,1), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(2,2), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(2,3), 1.0-2*z, 1e-14);
|
||||
EXPECT_NEAR(ten(3,1), 1.0-2*x, 1e-14);
|
||||
EXPECT_NEAR(ten(3,2), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(3,3), 0.0, 1e-14);
|
||||
|
||||
SymmTensor sten = f3D_sten(X);
|
||||
EXPECT_NEAR(sten(1,1), 0.0, 1e-14);
|
||||
EXPECT_NEAR(sten(1,2), 1.0-2*x, 1e-14);
|
||||
EXPECT_NEAR(sten(1,3), 1.0-2*y, 1e-14);
|
||||
EXPECT_NEAR(sten(2,1), 1.0-2*x, 1e-14);
|
||||
EXPECT_NEAR(sten(2,2), 0.0, 1e-14);
|
||||
EXPECT_NEAR(sten(2,3), 1.0-2*x, 1e-14);
|
||||
EXPECT_NEAR(sten(3,1), 1.0-2*y, 1e-14);
|
||||
EXPECT_NEAR(sten(3,2), 1.0-2*x, 1e-14);
|
||||
EXPECT_NEAR(sten(3,3), 0.0, 1e-14);
|
||||
}
|
||||
|
||||
|
||||
TEST(TestFieldFunctions, 3D2PLR)
|
||||
{
|
||||
FieldFunction f3D_scalar("src/Utility/Test/refdata/Field3D-2PLR",
|
||||
"Stokes-2", "v", 2, 0);
|
||||
|
||||
VecFieldFunction f3D_vec("src/Utility/Test/refdata/Field3D-2PLR",
|
||||
"Stokes-1", "u", 2, 0);
|
||||
|
||||
TensorFieldFunction f3D_ten("src/Utility/Test/refdata/Field3D-2PLR",
|
||||
"Stokes-1", "u_x,x|u_y,x|u_z,x|u_x,y|u_y,y|u_z,y|u_x,z|u_x,y|u_z,z", 2, 0);
|
||||
|
||||
STensorFieldFunction f3D_sten("src/Utility/Test/refdata/Field3D-2PLR",
|
||||
"Stokes-1", "u_x,x|u_y,y|u_z,z|u_x,y|u_y,z|u_z,x", 2, 0);
|
||||
|
||||
double param[3] = {0.25, 0.25, 0.25};
|
||||
Vec4 X(param);
|
||||
for (size_t pIdx = 0; pIdx < 2; ++pIdx) {
|
||||
f3D_scalar.initPatch(pIdx);
|
||||
double scal = f3D_scalar(X);
|
||||
double x = 0.25 + pIdx;
|
||||
double y = 0.5;
|
||||
double z = 0.5;
|
||||
EXPECT_NEAR(scal, 1.0-2*x, 1e-14);
|
||||
f3D_vec.initPatch(pIdx);
|
||||
Vec3 vec = f3D_vec(X);
|
||||
EXPECT_NEAR(vec[0], y*(1.0-y), 1e-14); // y*(1-y)
|
||||
EXPECT_NEAR(vec[1], z*(1.0-z), 1e-14); // z*(1-z)
|
||||
EXPECT_NEAR(vec[2], x*(1.0-x), 1e-14); // x*(1-x)
|
||||
f3D_ten.initPatch(pIdx);
|
||||
Tensor ten = f3D_ten(X);
|
||||
EXPECT_NEAR(ten(1,1), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(1,2), 1.0-2*y, 1e-14);
|
||||
EXPECT_NEAR(ten(1,3), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(2,1), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(2,2), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(2,3), 1.0-2*z, 1e-14);
|
||||
EXPECT_NEAR(ten(3,1), 1.0-2*x, 1e-14);
|
||||
EXPECT_NEAR(ten(3,2), 0.0, 1e-14);
|
||||
EXPECT_NEAR(ten(3,3), 0.0, 1e-14);
|
||||
|
||||
f3D_sten.initPatch(pIdx);
|
||||
SymmTensor sten = f3D_sten(X);
|
||||
EXPECT_NEAR(sten(1,1), 0.0, 1e-14);
|
||||
EXPECT_NEAR(sten(1,2), 1.0-2*y, 1e-14);
|
||||
EXPECT_NEAR(sten(1,3), 1.0-2*x, 1e-14);
|
||||
EXPECT_NEAR(sten(2,1), 1.0-2*y, 1e-14);
|
||||
EXPECT_NEAR(sten(2,2), 0.0, 1e-14);
|
||||
EXPECT_NEAR(sten(2,3), 1.0-2*z, 1e-14);
|
||||
EXPECT_NEAR(sten(3,1), 1.0-2*x, 1e-14);
|
||||
EXPECT_NEAR(sten(3,2), 1.0-2*z, 1e-14);
|
||||
EXPECT_NEAR(sten(3,3), 0.0, 1e-14);
|
||||
}
|
||||
}
|
||||
BIN
src/Utility/Test/refdata/Field2D-1P.hdf5
Normal file
BIN
src/Utility/Test/refdata/Field2D-1P.hdf5
Normal file
Binary file not shown.
10
src/Utility/Test/refdata/Field2D-1P.xml
Normal file
10
src/Utility/Test/refdata/Field2D-1P.xml
Normal file
@ -0,0 +1,10 @@
|
||||
<info>
|
||||
<entry name="u" description="primary" type="field" basis="Stokes-1" patches="1" components="2" />
|
||||
<entry name="v" description="primary" type="field" basis="Stokes-2" patches="1" components="1" />
|
||||
<entry name="u_x,x" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_y,x" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_x,y" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_y,y" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<levels>0</levels>
|
||||
<timestep constant="1" interval="1">1.000000</timestep>
|
||||
</info>
|
||||
BIN
src/Utility/Test/refdata/Field2D-1PLR.hdf5
Normal file
BIN
src/Utility/Test/refdata/Field2D-1PLR.hdf5
Normal file
Binary file not shown.
10
src/Utility/Test/refdata/Field2D-1PLR.xml
Normal file
10
src/Utility/Test/refdata/Field2D-1PLR.xml
Normal file
@ -0,0 +1,10 @@
|
||||
<info>
|
||||
<entry name="u" description="primary" type="field" basis="Stokes-1" patches="1" components="2" />
|
||||
<entry name="v" description="primary" type="field" basis="Stokes-2" patches="1" components="1" />
|
||||
<entry name="u_x,x" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_y,x" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_x,y" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_y,y" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<levels>0</levels>
|
||||
<timestep constant="1" interval="1">1.000000</timestep>
|
||||
</info>
|
||||
BIN
src/Utility/Test/refdata/Field2D-2P.hdf5
Normal file
BIN
src/Utility/Test/refdata/Field2D-2P.hdf5
Normal file
Binary file not shown.
10
src/Utility/Test/refdata/Field2D-2P.xml
Normal file
10
src/Utility/Test/refdata/Field2D-2P.xml
Normal file
@ -0,0 +1,10 @@
|
||||
<info>
|
||||
<entry name="u" description="primary" type="field" basis="Stokes-1" patches="2" components="2" />
|
||||
<entry name="v" description="primary" type="field" basis="Stokes-2" patches="2" components="1" />
|
||||
<entry name="u_x,x" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_y,x" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_x,y" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_y,y" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<levels>0</levels>
|
||||
<timestep constant="1" interval="1">1.000000</timestep>
|
||||
</info>
|
||||
BIN
src/Utility/Test/refdata/Field2D-2PLR.hdf5
Normal file
BIN
src/Utility/Test/refdata/Field2D-2PLR.hdf5
Normal file
Binary file not shown.
10
src/Utility/Test/refdata/Field2D-2PLR.xml
Normal file
10
src/Utility/Test/refdata/Field2D-2PLR.xml
Normal file
@ -0,0 +1,10 @@
|
||||
<info>
|
||||
<entry name="u" description="primary" type="field" basis="Stokes-1" patches="2" components="2" />
|
||||
<entry name="v" description="primary" type="field" basis="Stokes-2" patches="2" components="1" />
|
||||
<entry name="u_x,x" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_y,x" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_x,y" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_y,y" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<levels>0</levels>
|
||||
<timestep constant="1" interval="1">1.000000</timestep>
|
||||
</info>
|
||||
BIN
src/Utility/Test/refdata/Field3D-1P.hdf5
Normal file
BIN
src/Utility/Test/refdata/Field3D-1P.hdf5
Normal file
Binary file not shown.
15
src/Utility/Test/refdata/Field3D-1P.xml
Normal file
15
src/Utility/Test/refdata/Field3D-1P.xml
Normal file
@ -0,0 +1,15 @@
|
||||
<info>
|
||||
<entry name="u" description="primary" type="field" basis="Stokes-1" patches="1" components="3" />
|
||||
<entry name="v" description="primary" type="field" basis="Stokes-2" patches="1" components="1" />
|
||||
<entry name="u_x,x" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_y,x" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_z,x" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_x,y" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_y,y" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_z,y" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_x,z" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_y,z" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_z,z" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<levels>0</levels>
|
||||
<timestep constant="1" interval="1">1.000000</timestep>
|
||||
</info>
|
||||
BIN
src/Utility/Test/refdata/Field3D-1PLR.hdf5
Normal file
BIN
src/Utility/Test/refdata/Field3D-1PLR.hdf5
Normal file
Binary file not shown.
15
src/Utility/Test/refdata/Field3D-1PLR.xml
Normal file
15
src/Utility/Test/refdata/Field3D-1PLR.xml
Normal file
@ -0,0 +1,15 @@
|
||||
<info>
|
||||
<entry name="u" description="primary" type="field" basis="Stokes-1" patches="1" components="3" />
|
||||
<entry name="v" description="primary" type="field" basis="Stokes-2" patches="1" components="1" />
|
||||
<entry name="u_x,x" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_y,x" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_z,x" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_x,y" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_y,y" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_z,y" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_x,z" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_y,z" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<entry name="u_z,z" description="secondary" type="field" basis="Stokes-1" patches="1" components="1" />
|
||||
<levels>0</levels>
|
||||
<timestep constant="1" interval="1">1.000000</timestep>
|
||||
</info>
|
||||
BIN
src/Utility/Test/refdata/Field3D-2P.hdf5
Normal file
BIN
src/Utility/Test/refdata/Field3D-2P.hdf5
Normal file
Binary file not shown.
15
src/Utility/Test/refdata/Field3D-2P.xml
Normal file
15
src/Utility/Test/refdata/Field3D-2P.xml
Normal file
@ -0,0 +1,15 @@
|
||||
<info>
|
||||
<entry name="u" description="primary" type="field" basis="Stokes-1" patches="2" components="3" />
|
||||
<entry name="v" description="primary" type="field" basis="Stokes-2" patches="2" components="1" />
|
||||
<entry name="u_x,x" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_y,x" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_z,x" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_x,y" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_y,y" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_z,y" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_x,z" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_y,z" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_z,z" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<levels>0</levels>
|
||||
<timestep constant="1" interval="1">1.000000</timestep>
|
||||
</info>
|
||||
BIN
src/Utility/Test/refdata/Field3D-2PLR.hdf5
Normal file
BIN
src/Utility/Test/refdata/Field3D-2PLR.hdf5
Normal file
Binary file not shown.
15
src/Utility/Test/refdata/Field3D-2PLR.xml
Normal file
15
src/Utility/Test/refdata/Field3D-2PLR.xml
Normal file
@ -0,0 +1,15 @@
|
||||
<info>
|
||||
<entry name="u" description="primary" type="field" basis="Stokes-1" patches="2" components="3" />
|
||||
<entry name="v" description="primary" type="field" basis="Stokes-2" patches="2" components="1" />
|
||||
<entry name="u_x,x" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_y,x" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_z,x" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_x,y" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_y,y" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_z,y" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_x,z" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_y,z" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<entry name="u_z,z" description="secondary" type="field" basis="Stokes-1" patches="2" components="1" />
|
||||
<levels>0</levels>
|
||||
<timestep constant="1" interval="1">1.000000</timestep>
|
||||
</info>
|
||||
@ -192,22 +192,28 @@ class Vec4 : public Vec3
|
||||
public:
|
||||
Real t; //!< The time coordinate
|
||||
int idx; //!< Nodal point index
|
||||
const Real* u; //!< Spline parameters of point
|
||||
|
||||
//! \brief Default constructor creating a point at origin.
|
||||
Vec4() { t = 0.0; idx = -1; }
|
||||
Vec4(const Real* par = nullptr) { t = 0.0; idx = -1; u = par; }
|
||||
//! \brief Constructor creating a point at the specified location.
|
||||
Vec4(Real X, Real Y, Real Z, Real T = 0.0) : Vec3(X,Y,Z) { t = T; idx = -1; }
|
||||
Vec4(Real X, Real Y, Real Z, Real T = 0.0) : Vec3(X,Y,Z)
|
||||
{ t = T; idx = -1; u = nullptr; }
|
||||
//! \brief Constructor creating a point at the specified location.
|
||||
Vec4(const Vec3& X, Real T = 0.0, int id = -1) : Vec3(X) { t = T; idx = id; }
|
||||
Vec4(const Vec3& X, Real T = 0.0, int id = -1) : Vec3(X)
|
||||
{ t = T; idx = id; u = nullptr; }
|
||||
//! \brief Constructor creating a point at the specified location.
|
||||
Vec4(const Vec3& X, Real T, const Real* par) : Vec3(X)
|
||||
{ t = T; idx = -1; u = par; }
|
||||
//! \brief Constructor creating a point from the given \a std::vector.
|
||||
Vec4(const std::vector<Real>& X) : Vec3(X)
|
||||
{ t = X.size() > 3 ? X[3] : 0.0; idx = -1; }
|
||||
{ t = X.size() > 3 ? X[3] : 0.0; idx = -1; u = nullptr; }
|
||||
//! \brief Copy constructor.
|
||||
Vec4(const Vec4& X) : Vec3(X) { t = X.t; idx = X.idx; }
|
||||
Vec4(const Vec4& X) : Vec3(X) { t = X.t; idx = X.idx; u = X.u; }
|
||||
|
||||
//! \brief Assignment operator.
|
||||
Vec4& operator=(const Vec4& X)
|
||||
{ x = X.x; y = X.y; z = X.z; t = X.t; idx = X.idx; return *this; }
|
||||
{ x = X.x; y = X.y; z = X.z; t = X.t; idx = X.idx; u = X.u; return *this; }
|
||||
//! \brief Overloaded assignment operator.
|
||||
Vec4& operator=(Real val) { x = y = z = val; return *this; }
|
||||
|
||||
@ -216,7 +222,7 @@ public:
|
||||
{
|
||||
x = X.x; y = X.y; z = X.z;
|
||||
const Vec4* x4 = dynamic_cast<const Vec4*>(&X);
|
||||
if (x4) { t = x4->t; idx = x4->idx; }
|
||||
if (x4) { t = x4->t; idx = x4->idx; u = x4->u; }
|
||||
return *this;
|
||||
}
|
||||
|
||||
@ -225,6 +231,7 @@ public:
|
||||
{
|
||||
if (idx >= 0) os <<"(idx="<< idx <<") ";
|
||||
this->Vec3::print(os,tol);
|
||||
if (u) os <<" ("<< u[0] <<" "<< u[1] <<" "<< u[2] <<")";
|
||||
if (t > 0.0) os <<" "<< t;
|
||||
return os;
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user