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

Change the Lagrange::computeBasis interface such that derval may be omitted when not needed. Also alow for input coordinates extrapolating the Gauss points.

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1107 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo
2011-08-23 14:51:00 +00:00
committed by Knut Morten Okstad
parent 272ec9192e
commit bd17c61884
7 changed files with 128 additions and 118 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/ASM/ASMs1DSpec.C
View File

@@ -119,7 +119,7 @@ bool ASMs1DSpec::integrate (Integrand& integrand,
dNdu.fillColumn(1,D1.getRow(i+1));
}
else
if (!Lagrange::computeBasis(fe.N,dNdu,points1,xg1[i]))
if (!Lagrange::computeBasis(fe.N,&dNdu,points1,xg1[i]))
return false;
// Compute Jacobian inverse of coordinate mapping and derivatives

48
src/ASM/Lagrange.C
View File

@@ -16,15 +16,6 @@
bool Lagrange::evalPol (int polnum, double xi, double& retval) const
{
// Check if evaluation point is inside the range of interpolation points
if (xi < points.front() || xi > points.back())
{
std::cerr <<" *** Lagrange::evalPol: Evaluation point out of range: "
<< xi <<", must be in the interval ["<< points.front()
<<","<< points.back() <<"]"<< std::endl;
return false;
}
// Degree of polynomials = number of polynomials
size_t degree = points.size();
@@ -86,7 +77,22 @@ bool Lagrange::evalDer (int polnum, double xi, double& retval) const
}
bool Lagrange::computeBasis (Vector& val,
bool Lagrange::computeBasis (RealArray& val,
int p1, double x1,
int p2, double x2,
int p3, double x3)
{
// Define the Lagrangian interpolation points
RealArray points1(p1), points2(p2), points3(p3);
for (int i = 0; i < p1; i++) points1[i] = -1.0 + (i+i)/double(p1-1);
for (int j = 0; j < p2; j++) points2[j] = -1.0 + (j+j)/double(p2-1);
for (int k = 0; k < p3; k++) points3[k] = -1.0 + (k+k)/double(p3-1);
return Lagrange::computeBasis(val,0,points1,x1,points2,x2,points3,x3);
}
bool Lagrange::computeBasis (RealArray& val,
Matrix& derval,
int p1, double x1,
int p2, double x2,
@@ -98,12 +104,12 @@ bool Lagrange::computeBasis (Vector& val,
for (int j = 0; j < p2; j++) points2[j] = -1.0 + (j+j)/double(p2-1);
for (int k = 0; k < p3; k++) points3[k] = -1.0 + (k+k)/double(p3-1);
return Lagrange::computeBasis(val,derval,points1,x1,points2,x2,points3,x3);
return Lagrange::computeBasis(val,&derval,points1,x1,points2,x2,points3,x3);
}
bool Lagrange::computeBasis (Vector& val,
Matrix& derval,
bool Lagrange::computeBasis (RealArray& val,
Matrix* derval,
const RealArray& points1, double x1,
const RealArray& points2, double x2,
const RealArray& points3, double x3)
@@ -144,11 +150,15 @@ bool Lagrange::computeBasis (Vector& val,
// Evaluate values of the 3-dimensional basis functions in the point x
val.resize(nen);
size_t count = 1;
size_t count = 0;
for (k = 0; k < n3; k++)
for (j = 0; j < n2; j++)
for (i = 0; i < n1; i++, count++)
val(count) = tempval1[i]*tempval2[j]*tempval3[k];
val[count] = tempval1[i]*tempval2[j]*tempval3[k];
if (!derval) return true;
Matrix& deriv = *derval;
// Vectors of derivative values for the one dimensional polynomials
RealArray tempder1(p1), tempder2(p2), tempder3(p3);
@@ -167,16 +177,16 @@ bool Lagrange::computeBasis (Vector& val,
return false;
// Evaluate derivatives of the 3-dimensional basis functions in the point x
derval.resize(nen, p3 > 0 ? 3 : (p2 > 0 ? 2 : 1));
derval->resize(nen, p3 > 0 ? 3 : (p2 > 0 ? 2 : 1));
count = 1;
for (k = 0; k < n3; k++)
for (j = 0; j < n2; j++)
for (i = 0; i < n1; i++, count++)
{
if (p1 > 0) derval(count,1) = tempder1[i]*tempval2[j]*tempval3[k];
if (p2 > 0) derval(count,2) = tempval1[i]*tempder2[j]*tempval3[k];
if (p3 > 0) derval(count,3) = tempval1[i]*tempval2[j]*tempder3[k];
if (p1 > 0) deriv(count,1) = tempder1[i]*tempval2[j]*tempval3[k];
if (p2 > 0) deriv(count,2) = tempval1[i]*tempder2[j]*tempval3[k];
if (p3 > 0) deriv(count,3) = tempval1[i]*tempval2[j]*tempder3[k];
}
return true;

23
src/ASM/Lagrange.h
View File

@@ -41,7 +41,6 @@ public:
//! \brief Evaluates a 1D, 2D or 3D Lagrangian basis at a given point.
//! \param[out] val Values of all basis functions
//! \param[out] derval Derivatives of all basis functions
//! \param[in] p1 Polynomial degree in first parameter direction
//! \param[in] x1 Natural coordinate in first parameter direction
//! \param[in] p2 Polynomial degree in second parameter direction
@@ -53,14 +52,30 @@ public:
//! If \a p2 is zero, a 1D basis is assumed. Otherwise,
//! if \a p3 is zero, a 2D basis is assumed. If \a p1, \a p2, and \a p3
//! all are non-zero, a 3D basis is assumed.
static bool computeBasis (Vector& val, Matrix& derval,
int p1 = 0, double x1 = 0.0,
static bool computeBasis (RealArray& val, int p1, double x1,
int p2 = 0, double x2 = 0.0,
int p3 = 0, double x3 = 0.0);
//! \brief Evaluates a 1D, 2D or 3D Lagrangian basis at a given point.
//! \param[out] val Values of all basis functions
//! \param[out] derval Derivatives of all basis functions
//! \param[in] p1 Polynomial degree in first parameter direction
//! \param[in] x1 Natural coordinate in first parameter direction
//! \param[in] p2 Polynomial degree in second parameter direction
//! \param[in] x2 Natural coordinate in second parameter direction
//! \param[in] p3 Polynomial degree in third parameter direction
//! \param[in] x3 Natural coordinate in third parameter direction
//!
//! \details If \a p2 is zero, a 1D basis is assumed. Otherwise,
//! if \a p3 is zero, a 2D basis is assumed. If \a p1, \a p2, and \a p3
//! all are non-zero, a 3D basis is assumed.
static bool computeBasis (RealArray& val, Matrix& derval, int p1, double x1,
int p2 = 0, double x2 = 0.0,
int p3 = 0, double x3 = 0.0);
//! \brief Evaluates a 1D, 2D or 3D Lagrangian basis at a given point.
//! \param[out] val Values of all basis functions
//! \param[out] derval Pointer to derivatives of all basis functions
//! \param[in] p1 Natural point coordinates in first parameter direction
//! \param[in] x1 Natural coordinate in first parameter direction
//! \param[in] p2 Natural point coordinates in second parameter direction
@@ -72,7 +87,7 @@ public:
//! If \a p2 is empty, a 1D basis is assumed. Otherwise,
//! if \a p3 is empty, a 2D basis is assumed. If \a p1, \a p2, and \a p3
//! all are non-empty, a 3D basis is assumed.
static bool computeBasis (Vector& val, Matrix& derval,
static bool computeBasis (RealArray& val, Matrix* derval,
const RealArray& p1, double x1,
const RealArray& p2 = RealArray(), double x2 = 0.0,
const RealArray& p3 = RealArray(), double x3 = 0.0);

37
src/ASM/LagrangeField2D.C
View File

@@ -35,22 +35,22 @@ double LagrangeField2D::valueNode(int node) const
double LagrangeField2D::valueFE(const FiniteElement& fe) const
{
Vector N;
Matrix dNdu;
Lagrange::computeBasis(N,dNdu,p1,fe.xi,p2,fe.eta);
Lagrange::computeBasis(N,p1,fe.xi,p2,fe.eta);
const int nel1 = (n1-1)/p1;
div_t divresult = div(fe.iel,nel1);
const int iel1 = divresult.rem;
const int iel2 = divresult.quot;
int iel1 = divresult.rem;
int iel2 = divresult.quot;
const int node1 = p1*iel1-1;
const int node2 = p2*iel2-1;
int node;
int node, i, j;
int locNode = 1;
double value = 0.0;
for (int j = node2;j <= node2+p2;j++)
for (int i = node1;i <= node1+p1;i++, locNode++) {
for (j = node2; j <= node2+p2; j++)
for (i = node1; i <= node1+p1; i++, locNode++)
{
node = (j-1)*n1 + i;
value += values(node)*N(locNode);
}
@@ -78,18 +78,19 @@ bool LagrangeField2D::gradFE(const FiniteElement& fe, Vector& grad) const
const int nel1 = (n1-1)/p1;
div_t divresult = div(fe.iel,nel1);
const int iel1 = divresult.rem;
const int iel2 = divresult.quot;
int iel1 = divresult.rem;
int iel2 = divresult.quot;
const int node1 = p1*iel1-1;
const int node2 = p2*iel2-1;
const int nen = (p1+1)*(p2+1);
Matrix Xnod(nsd,nen);
int node;
int node, i, j;
int locNode = 1;
for (int j = node2;j <= node2+p2;j++)
for (int i = node1;i <= node1+p1;i++, locNode++) {
for (j = node2; j <= node2+p2; j++)
for (i = node1; i <= node1+p1; i++, locNode++)
{
node = (j-1)*n1 + i;
Xnod.fillColumn(locNode,coord.getColumn(node));
}
@@ -97,14 +98,12 @@ bool LagrangeField2D::gradFE(const FiniteElement& fe, Vector& grad) const
Matrix Jac, dNdX;
utl::Jacobian(Jac,dNdX,Xnod,dNdu,false);
double value;
locNode = 1;
for (int j = node2;j <= node2+p2;j++)
for (int i = node1;i <= node1+p1;i++, locNode++) {
for (j = node2; j <= node2+p2; j++)
for (i = node1; i <= node1+p1; i++, locNode++)
{
node = (j-1)*n1 + i;
value = values(node);
for (int k = 1;k <= nsd;k++)
grad(k) += value*dNdX(locNode,k);
grad.add(dNdX.getColumn(locNode),values(node));
}
return true;

30
src/ASM/LagrangeField3D.C
View File

@@ -36,8 +36,7 @@ double LagrangeField3D::valueNode(int node) const
double LagrangeField3D::valueFE(const FiniteElement& fe) const
{
Vector N;
Matrix dNdu;
Lagrange::computeBasis(N,dNdu,p1,fe.xi,p2,fe.eta,p3,fe.zeta);
Lagrange::computeBasis(N,p1,fe.xi,p2,fe.eta,p3,fe.zeta);
const int nel1 = (n1-1)/p1;
const int nel2 = (n2-1)/p2;
@@ -52,16 +51,16 @@ double LagrangeField3D::valueFE(const FiniteElement& fe) const
const int node2 = p2*iel2-1;
const int node3 = p3*iel3-1;
int node;
int node, i, j, k;
int locNode = 1;
double value = 0.0;
for (int k = node3;k <= node3+p3;k++)
for (int j = node2;j <= node2+p2;j++)
for (int i = node1;i <= node1+p1;i++, locNode++) {
for (k = node3; k <= node3+p3; k++)
for (j = node2; j <= node2+p2; j++)
for (i = node1; i <= node1+p1; i++, locNode++)
{
node = (k-1)*n1*n2 + (j-1)*n1 + i;
value += values(node)*N(locNode);
}
return value;
}
@@ -95,15 +94,15 @@ bool LagrangeField3D::gradFE(const FiniteElement& fe, Vector& grad) const
const int node1 = p1*iel1-1;
const int node2 = p2*iel2-1;
const int node3 = p3*iel3-1;
const int nen = (p1+1)*(p2+1)*(p3+1);
Matrix Xnod(nsd,nen);
int node, i, j, k;
int locNode = 1;
for (k = node3;k <= node3+p3;k++)
for (j = node2;j <= node2+p2;j++)
for (i = node1;i <= node1+p1;i++, locNode++) {
for (k = node3; k <= node3+p3; k++)
for (j = node2; j <= node2+p2; j++)
for (i = node1; i <= node1+p1; i++, locNode++) {
node = (k-1)*n1*n2 + (j-1)*n1 + i;
Xnod.fillColumn(locNode,coord.getColumn(node));
}
@@ -111,15 +110,12 @@ bool LagrangeField3D::gradFE(const FiniteElement& fe, Vector& grad) const
Matrix Jac, dNdX;
utl::Jacobian(Jac,dNdX,Xnod,dNdu,false);
double value;
locNode = 1;
for (k = node3;k <= node3+p3;k++)
for (j = node2;j <= node2+p2;j++)
for (k = node3; k <= node3+p3; k++)
for (j = node2; j <= node2+p2; j++)
for (i = node1;i <= node1+p1;i++, locNode++) {
node = (k-1)*n1*n2 + (j-1)*n1 + i;
value = values(node);
for (int k = 1;k <= nsd;k++)
grad(k) += value*dNdX(locNode,k);
grad.add(dNdX.getColumn(locNode),values(node));
}
return true;

53
src/ASM/LagrangeFields2D.C
View File

@@ -41,30 +41,28 @@ bool LagrangeFields2D::valueNode(int node, Vector& vals) const
bool LagrangeFields2D::valueFE(const FiniteElement& fe, Vector& vals) const
{
vals.resize(nf,0.0);
Vector N;
Matrix dNdu;
Lagrange::computeBasis(N,dNdu,p1,fe.xi,p2,fe.eta);
Lagrange::computeBasis(N,p1,fe.xi,p2,fe.eta);
const int nel1 = (n1-1)/p1;
div_t divresult = div(fe.iel,nel1);
const int iel1 = divresult.rem;
const int iel2 = divresult.quot;
int iel1 = divresult.rem;
int iel2 = divresult.quot;
const int node1 = p1*iel1-1;
const int node2 = p2*iel2-1;
int node, dof;
int i, j, k, dof;
int locNode = 1;
double value;
for (int j = node2;j <= node2+p2;j++)
for (int i = node1;i <= node1+p1;i++, locNode++) {
node = (j-1)*n1 + i;
for (j = node2; j <= node2+p2; j++)
for (i = node1; i <= node1+p1; i++, locNode++)
{
dof = nf*((j-1)*n1 + i-1) + 1;
value = N(locNode);
for (int k = 1;k <= nf;k++) {
dof = nf*(node-1) + k;
for (k = 1; k <= nf; k++, dof++)
vals(k) += values(dof)*value;
}
}
return true;
@@ -91,18 +89,19 @@ bool LagrangeFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const
const int nel1 = (n1-1)/p1;
div_t divresult = div(fe.iel,nel1);
const int iel1 = divresult.rem;
const int iel2 = divresult.quot;
int iel1 = divresult.rem;
int iel2 = divresult.quot;
const int node1 = p1*iel1-1;
const int node2 = p2*iel2-1;
const int nen = (p1+1)*(p2+1);
Matrix Xnod(nsd,nen);
int node;
int node, dof, i, j;
int locNode = 1;
for (int j = node2;j <= node2+p2;j++)
for (int i = node1;i <= node1+p1;i++, locNode++) {
for (j = node2; j <= node2+p2; j++)
for (i = node1; i <= node1+p1; i++, locNode++)
{
node = (j-1)*n1 + i;
Xnod.fillColumn(locNode,coord.getColumn(node));
}
@@ -110,18 +109,14 @@ bool LagrangeFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const
Matrix Jac, dNdX;
utl::Jacobian(Jac,dNdX,Xnod,dNdu,false);
double value;
int dof;
locNode = 1;
for (int j = node2;j <= node2+p2;j++)
for (int i = node1;i <= node1+p1;i++, locNode++) {
node = (j-1)*n1 + i;
for (int k = 1;k <= nf;k++) {
dof = nf*(node-1) + k;
value = values(dof);
for (int l = 1;l <= nsd;l++)
grad(k,l) += value*dNdX(locNode,l);
}
for (j = node2; j <= node2+p2; j++)
for (i = node1; i <= node1+p1; i++, locNode++)
{
dof = nf*((j-1)*n1 + i-1) + 1;
for (int k = 1; k <= nf; k++, dof++)
for (int l = 1; l <= nsd; l++)
grad(k,l) += values(dof)*dNdX(locNode,l);
}
return true;

53
src/ASM/LagrangeFields3D.C
View File

@@ -43,10 +43,9 @@ bool LagrangeFields3D::valueNode(int node, Vector& vals) const
bool LagrangeFields3D::valueFE(const FiniteElement& fe, Vector& vals) const
{
vals.resize(nf,0.0);
Vector N;
Matrix dNdu;
Lagrange::computeBasis(N,dNdu,p1,fe.xi,p2,fe.eta,p3,fe.zeta);
Lagrange::computeBasis(N,p1,fe.xi,p2,fe.eta,p3,fe.zeta);
const int nel1 = (n1-1)/p1;
const int nel2 = (n2-1)/p2;
@@ -61,20 +60,19 @@ bool LagrangeFields3D::valueFE(const FiniteElement& fe, Vector& vals) const
const int node2 = p2*iel2-1;
const int node3 = p3*iel3-1;
int node, dof;
int i, j, k, l, dof;
int locNode = 1;
double value;
for (int k = node3;k <= node3+p3;k++)
for (int j = node2;j <= node2+p2;j++)
for (int i = node1;i <= node1+p1;i++, locNode++) {
node = (k-1)*n1*n2 + (j-1)*n1 + i;
for (k = node3; k <= node3+p3; k++)
for (j = node2; j <= node2+p2; j++)
for (i = node1; i <= node1+p1; i++, locNode++)
{
dof = nf*(n1*(n2*(k-1) + j-1) + i-1) + 1;
value = N(locNode);
for (int l = 1;l <= nf;l++) {
dof = nf*(node-1)+l;
for (l = 1; l <= nf; l++, dof++)
vals(l) += values(dof)*value;
}
}
return true;
}
@@ -108,15 +106,16 @@ bool LagrangeFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const
const int node1 = p1*iel1-1;
const int node2 = p2*iel2-1;
const int node3 = p3*iel3-1;
const int nen = (p1+1)*(p2+1)*(p3+1);
Matrix Xnod(nsd,nen);
int node;
int node, dof, i, j, k;
int locNode = 1;
for (int k = node3;k <= node3+p3;k++)
for (int j = node2;j <= node2+p2;j++)
for (int i = node1;i <= node1+p1;i++, locNode++) {
for (k = node3; k <= node3+p3; k++)
for (j = node2; j <= node2+p2; j++)
for (i = node1; i <= node1+p1; i++, locNode++)
{
node = (k-1)*n1*n2 + (j-1)*n1 + i;
Xnod.fillColumn(locNode,coord.getColumn(node));
}
@@ -124,19 +123,15 @@ bool LagrangeFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const
Matrix Jac, dNdX;
utl::Jacobian(Jac,dNdX,Xnod,dNdu,false);
int dof;
double value;
locNode = 1;
for (int k = node3;k <= node3+p3;k++)
for (int j = node2;j <= node2+p2;j++)
for (int i = node1;i <= node1+p1;i++, locNode++) {
node = (k-1)*n1*n2 + (j-1)*n1 + i;
for (int m = 1;m <= nf;m++) {
dof = nf*(node-1) + m;
value = values(dof);
for (int n = 1;n <= nsd;n++)
grad(m,n) += value*dNdX(locNode,n);
}
for (k = node3; k <= node3+p3; k++)
for (j = node2; j <= node2+p2; j++)
for (i = node1; i <= node1+p1; i++, locNode++)
{
dof = nf*(n1*(n2*(k-1) + j-1) + i-1) + 1;
for (int m = 1; m <= nf; m++, dof++)
for (int n = 1; n <= nsd; n++)
grad(m,n) += values(dof)*dNdX(locNode,n);
}
return true;