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

Field definition for Lagrange finite element fields

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@1051 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
rho 2011-06-10 13:31:18 +00:00 committed by Knut Morten Okstad
parent d451094a42
commit aae2a49143
8 changed files with 819 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

118
src/ASM/LagrangeField2D.C Normal file
View File

@ -0,0 +1,118 @@
//==============================================================================
//!
//! \file LagrangeField2D.C
//!
//! \date Jun 8 2011
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Class for Lagrange-based finite element scalar field in 2D.
//!
//==============================================================================
#include "LagrangeField2D.h"
#include "FiniteElement.h"
#include "Lagrange.h"
#include "CoordinateMapping.h"
#include "Vec3.h"
LagrangeField2D::LagrangeField2D(Matrix X, int nx, int ny,
int px, int py, char* name)
: Field(2,name), coord(X), n1(nx), n2(ny), p1(px), p2(py)
{
nno = n1*n2;
nelm = (n1-1)*(n2-1)/(p1*p2);
}
double LagrangeField2D::valueNode(int node) const
{
return values(node);
}
double LagrangeField2D::valueFE(const FiniteElement& fe) const
{
Vector N;
Matrix dNdu;
Lagrange::computeBasis(N,dNdu,p1,fe.u,p2,fe.v);
const int nel1 = (n1-1)/p1;
div_t divresult = div(fe.iel,nel1);
const int iel1 = divresult.rem;
const int iel2 = divresult.quot;
const int node1 = p1*iel1-1;
const int node2 = p2*iel2-1;
int node;
int locNode = 1;
double value = 0.0;
for (int j = node2;j <= node2+p2;j++)
for (int i = node1;i <= node1+p1;i++, locNode++) {
node = (j-1)*n1 + i;
value += values(node)*N(locNode);
}
return value;
}
double LagrangeField2D::valueCoor(const Vec3& x) const
{
// Not implemented yet
return 0.0;
}
bool LagrangeField2D::gradFE(const FiniteElement& fe, Vector& grad) const
{
grad.resize(nsd,0.0);
Vector N;
Matrix dNdu;
if (!Lagrange::computeBasis(N,dNdu,p1,fe.u,p2,fe.v))
return false;
const int nel1 = (n1-1)/p1;
div_t divresult = div(fe.iel,nel1);
const int iel1 = divresult.rem;
const 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 locNode = 1;
for (int j = node2;j <= node2+p2;j++)
for (int i = node1;i <= node1+p1;i++, locNode++) {
node = (j-1)*n1 + i;
Xnod.fillColumn(locNode,coord.getColumn(node));
}
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++) {
node = (j-1)*n1 + i;
value = values(node);
for (int k = 1;k <= nsd;k++)
grad(k) += value*dNdX(locNode,k);
}
return true;
}
bool LagrangeField2D::gradCoor(const Vec3& x, Vector& grad) const
{
// Not implemented yet
return false;
}

71
src/ASM/LagrangeField2D.h Normal file
View File

@ -0,0 +1,71 @@
//==============================================================================
//!
//! \file LagrangeField2D.h
//!
//! \date June 8 2011
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Class for Lagrange-based finite element scalar fields in 2D.
//!
//==============================================================================
#ifndef _LAGRANGE_FIELD_2D_H
#define _LAGRANGE_FIELD_2D_H
#include "Field.h"
/*!
\brief Class for Lagrange-based finite element scalar fields in 2D.
\details This class implements the functions required to evaluate a 2D
Lagrange scalar field at a given point in parametrical or physical coordinates.
*/
class LagrangeField2D : public Field
{
public:
//! \brief The constructor sets the number of space dimensions and fields.
//! \param[in] X Matrix of nodel coordinates
//! \param[in] n1 Number of nodes in first parameter direction
//! \param[in] n2 Number of nodes in second parameter direction
//! \param[in] p1 Element order in first parameter direction
//! \param[in] p2 Element order in second parameter direction
//! \param[in] name Name of spline field
LagrangeField2D(Matrix X, int nx, int ny, int px, int py, char* name = NULL);
//! \brief Empty destructor.
virtual ~LagrangeField2D() {}
// Methods to compute field values
//================================
//! \brief Computes the value in a given node/control point.
//! \param[in] node Node number
double valueNode(int node) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
double valueFE(const FiniteElement& fe) const;
//! \brief Computed the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point
double valueCoor(const Vec3& x) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[out] grad Gradient of solution in a given local coordinate
bool gradFE(const FiniteElement& fe, Vector& grad) const;
//! \brief Computes the gradient for a given global/physical coordinate.
//! \param[in] x Global coordinate
//! \param[out] grad Gradient of solution in a given global coordinate
bool gradCoor(const Vec3& x, Vector& grad) const;
protected:
Matrix coord; //!< Matrix of nodel coordinates
int n1, n2; //!< Number of nodes in each parameter direction
int p1, p2; //!< Element order in each parameter direction
};
#endif

133
src/ASM/LagrangeField3D.C Normal file
View File

@ -0,0 +1,133 @@
//==============================================================================
//!
//! \file LagrangeField3D.C
//!
//! \date Jun 8 2011
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Class for Lagrange-based finite element scalar field in 3D.
//!
//==============================================================================
#include "LagrangeField3D.h"
#include "FiniteElement.h"
#include "Lagrange.h"
#include "CoordinateMapping.h"
#include "Vec3.h"
LagrangeField3D::LagrangeField3D(Matrix X, int nx, int ny, int nz,
int px, int py, int pz, char* name)
: Field(2,name), coord(X), n1(nx), n2(ny), n3(nz),
p1(px), p2(py), p3(pz)
{
nno = n1*n2*n3;
nelm = (n1-1)*(n2-1)*(n3-1)/(p1*p2*p3);
}
double LagrangeField3D::valueNode(int node) const
{
return values(node);
}
double LagrangeField3D::valueFE(const FiniteElement& fe) const
{
Vector N;
Matrix dNdu;
Lagrange::computeBasis(N,dNdu,p1,fe.u,p2,fe.v,p3,fe.w);
const int nel1 = (n1-1)/p1;
const int nel2 = (n2-1)/p2;
div_t divresult = div(fe.iel,nel1*nel2);
int iel2 = divresult.rem;
int iel3 = divresult.quot;
divresult = div(iel2,nel1);
int iel1 = divresult.rem;
iel2 = divresult.quot;
const int node1 = p1*iel1-1;
const int node2 = p2*iel2-1;
const int node3 = p3*iel3-1;
int node;
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++) {
node = (k-1)*n1*n2 + (j-1)*n1 + i;
value += values(node)*N(locNode);
}
return value;
}
double LagrangeField3D::valueCoor(const Vec3& x) const
{
// Not implemented yet
return 0.0;
}
bool LagrangeField3D::gradFE(const FiniteElement& fe, Vector& grad) const
{
grad.resize(nsd,0.0);
Vector N;
Matrix dNdu;
if (!Lagrange::computeBasis(N,dNdu,p1,fe.u,p2,fe.v,p3,fe.w))
return false;
const int nel1 = (n1-1)/p1;
const int nel2 = (n2-1)/p2;
div_t divresult = div(fe.iel,nel1*nel2);
int iel2 = divresult.rem;
int iel3 = divresult.quot;
divresult = div(iel2,nel1);
int iel1 = divresult.rem;
iel2 = divresult.quot;
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++) {
node = (k-1)*n1*n2 + (j-1)*n1 + i;
Xnod.fillColumn(locNode,coord.getColumn(node));
}
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 (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);
}
return true;
}
bool LagrangeField3D::gradCoor(const Vec3& x, Vector& grad) const
{
// Not implemented yet
return false;
}

68
src/ASM/LagrangeField3D.h Normal file
View File

@ -0,0 +1,68 @@
//==============================================================================
//!
//! \file LagrangeField3D.h
//!
//! \date Jun 8 2011
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Class for Lagrange-based finite element scalar fields in 3D.
//!
//==============================================================================
#ifndef _LAGRANGE_FIELD_3D_H
#define _LAGRANGE_FIELD_3D_H
#include "Field.h"
/*!
\brief Class for Lagrange-based finite element scalar fields in 3D.
\details This class implements the functions required to evaluate a 3D
Lagrange scalar field at a given point in parametrical or physical coordinates.
*/
class LagrangeField3D : public Field
{
public:
//! \brief The constructor sets the number of space dimensions and fields.
//! \param[in] geometry Spline volume geometry
//! \param[in] name Name of spline field
LagrangeField3D(Matrix X, int nx, int ny, int nz,
int px, int py, int pz, char* name = NULL);
//! \brief Empty destructor.
virtual ~LagrangeField3D() {}
// Methods to compute field values
//================================
//! \brief Computes the value in a given node/control point.
//! \param[in] node Node number
double valueNode(int node) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
double valueFE(const FiniteElement& fe) const;
//! \brief Computed the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point
double valueCoor(const Vec3& x) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[out] grad Gradient of solution in a given local coordinate
bool gradFE(const FiniteElement& fe, Vector& grad) const;
//! \brief Computes the gradient for a given global/physical coordinate.
//! \param[in] x Global coordinate
//! \param[out] grad Gradient of solution in a given global coordinate
bool gradCoor(const Vec3& x, Vector& grad) const;
protected:
Matrix coord; //!< Matrix of nodel coordinates
int n1, n2, n3; //!< Number of nodes in each parameter direction
int p1, p2, p3; //!< Element order in each parameter direction
};
#endif

135
src/ASM/LagrangeFields2D.C Normal file
View File

@ -0,0 +1,135 @@
//==============================================================================
//!
//! \file LagrangeFields2D.C
//!
//! \date Jun 8 2011
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Class for Lagrange-based finite element vector fields in 2D.
//!
//==============================================================================
#include "LagrangeFields2D.h"
#include "FiniteElement.h"
#include "Lagrange.h"
#include "CoordinateMapping.h"
#include "Vec3.h"
LagrangeFields2D::LagrangeFields2D(Matrix X, int nx, int ny,
int px, int py, char* name)
: Fields(2,name), coord(X), n1(nx), n2(ny), p1(px), p2(py)
{
nno = n1*n2;
nelm = (n1-1)*(n2-1)/(p1*p2);
// Number of fields set in fill
nf = 0;
}
bool LagrangeFields2D::valueNode(int node, Vector& vals) const
{
vals.resize(nf,0.0);
for (int i = 1;i <= nf;i++)
vals(i) = values(nf*(node-1)+i);
return true;
}
bool LagrangeFields2D::valueFE(const FiniteElement& fe, Vector& vals) const
{
vals.resize(nf,0.0);
Vector N;
Matrix dNdu;
Lagrange::computeBasis(N,dNdu,p1,fe.u,p2,fe.v);
const int nel1 = (n1-1)/p1;
div_t divresult = div(fe.iel,nel1);
const int iel1 = divresult.rem;
const int iel2 = divresult.quot;
const int node1 = p1*iel1-1;
const int node2 = p2*iel2-1;
int node, 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;
value = N(locNode);
for (int k = 1;k <= nf;k++) {
dof = nf*(node-1) + k;
vals(k) += values(dof)*value;
}
}
return true;
}
bool LagrangeFields2D::valueCoor(const Vec3& x, Vector& vals) const
{
// Not implemented yet
return false;
}
bool LagrangeFields2D::gradFE(const FiniteElement& fe, Matrix& grad) const
{
grad.resize(nf,nsd);
grad.fill(0.0);
Vector N;
Matrix dNdu;
if (!Lagrange::computeBasis(N,dNdu,p1,fe.u,p2,fe.v))
return false;
const int nel1 = (n1-1)/p1;
div_t divresult = div(fe.iel,nel1);
const int iel1 = divresult.rem;
const 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 locNode = 1;
for (int j = node2;j <= node2+p2;j++)
for (int i = node1;i <= node1+p1;i++, locNode++) {
node = (j-1)*n1 + i;
Xnod.fillColumn(locNode,coord.getColumn(node));
}
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);
}
}
return true;
}
bool LagrangeFields2D::gradCoor(const Vec3& x, Matrix& grad) const
{
// Not implemented yet
return false;
}

74
src/ASM/LagrangeFields2D.h Normal file
View File

@ -0,0 +1,74 @@
//==============================================================================
//!
//! \file LagrangeFields2D.h
//!
//! \date Jun 8 2011
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Class for Lagrange-based finite element vector fields in 2D.
//!
//==============================================================================
#ifndef _LAGRANGE_FIELDS_2D_H
#define _LAGRANGE_FIELDS_2D_H
#include "Fields.h"
/*!
\brief Class for Lagrange-based finite element vector fields in 2D.
\details This class implements the functions required to evaluate a 2D
Lagrange vector field at a given point in parametrical or physical coordinates.
*/
class LagrangeFields2D : public Fields
{
public:
//! \brief The constructor sets the field name.
//! \param[in] X Matrix of nodel coordinates
//! \param[in] n1 Number of nodes in first parameter direction
//! \param[in] n2 Number of nodes in second parameter direction
//! \param[in] p1 Element order in first parameter direction
//! \param[in] p2 Element order in second parameter direction
//! \param[in] name Name of spline field
LagrangeFields2D(Matrix X, int nx, int ny, int px, int py, char* name = NULL);
//! \brief Empty destructor.
virtual ~LagrangeFields2D() {}
// Methods to compute field values
//================================
//! \brief Computes the value in a given node/control point.
//! \param[in] node Node number
//! \param[out] vals Node values
bool valueNode(int node, Vector& vals) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[out] vals Values in local point in given element
bool valueFE(const FiniteElement& fe, Vector& vals) const;
//! \brief Computed the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point
//! \param[in] vals Values in given physical coordinate
bool valueCoor(const Vec3& x, Vector& vals) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[out] grad Gradient of solution in a given local coordinate
bool gradFE(const FiniteElement& fe, Matrix& grad) const;
//! \brief Computes the gradient for a given global/physical coordinate.
//! \param[in] x Global coordinate
//! \param[out] grad Gradient of solution in a given global coordinate
bool gradCoor(const Vec3& x, Matrix& grad) const;
protected:
Matrix coord; //!< Matrix of nodel coordinates
int n1, n2; //!< Number of nodes in each parameter direction
int p1, p2; //!< Element order in each parameter direction
};
#endif

150
src/ASM/LagrangeFields3D.C Normal file
View File

@ -0,0 +1,150 @@
//==============================================================================
//!
//! \file LagrangeFields3D.C
//!
//! \date Jun 8 2011
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Class for Lagrange-based finite element vector fields in 3D.
//!
//==============================================================================
#include "LagrangeFields3D.h"
#include "FiniteElement.h"
#include "Lagrange.h"
#include "CoordinateMapping.h"
#include "Vec3.h"
LagrangeFields3D::LagrangeFields3D(Matrix X, int nx, int ny, int nz,
int px, int py, int pz, char* name)
: Fields(3,name), coord(X), n1(nx), n2(ny), n3(nz),
p1(px), p2(py), p3(pz)
{
nno = n1*n2*n3;
nelm = (n1-1)*(n2-1)*(n3-1)/(p1*p2*p3);
// Number of fields set in fill
nf = 0;
}
bool LagrangeFields3D::valueNode(int node, Vector& vals) const
{
vals.resize(nf,0.0);
for (int i = 1;i <= nf;i++)
vals(i) = values(nf*(node-1)+i);
return true;
}
bool LagrangeFields3D::valueFE(const FiniteElement& fe, Vector& vals) const
{
vals.resize(nf,0.0);
Vector N;
Matrix dNdu;
Lagrange::computeBasis(N,dNdu,p1,fe.u,p2,fe.v,p3,fe.w);
const int nel1 = (n1-1)/p1;
const int nel2 = (n2-1)/p2;
div_t divresult = div(fe.iel,nel1*nel2);
int iel2 = divresult.rem;
int iel3 = divresult.quot;
divresult = div(iel2,nel1);
int iel1 = divresult.rem;
iel2 = divresult.quot;
const int node1 = p1*iel1-1;
const int node2 = p2*iel2-1;
const int node3 = p3*iel3-1;
int node, 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;
value = N(locNode);
for (int l = 1;l <= nf;l++) {
dof = nf*(node-1)+l;
vals(l) += values(dof)*value;
}
}
return true;
}
bool LagrangeFields3D::valueCoor(const Vec3& x, Vector& vals) const
{
// Not implemented yet
return false;
}
bool LagrangeFields3D::gradFE(const FiniteElement& fe, Matrix& grad) const
{
grad.resize(nf,nsd);
grad.fill(0.0);
Vector N;
Matrix dNdu;
if (!Lagrange::computeBasis(N,dNdu,p1,fe.u,p2,fe.v,p3,fe.w))
return false;
const int nel1 = (n1-1)/p1;
const int nel2 = (n2-1)/p2;
div_t divresult = div(fe.iel,nel1*nel2);
int iel2 = divresult.rem;
int iel3 = divresult.quot;
divresult = div(iel2,nel1);
int iel1 = divresult.rem;
iel2 = divresult.quot;
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 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;
Xnod.fillColumn(locNode,coord.getColumn(node));
}
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);
}
}
return true;
}
bool LagrangeFields3D::gradCoor(const Vec3& x, Matrix& grad) const
{
// Not implemented yet
return false;
}

70
src/ASM/LagrangeFields3D.h Normal file
View File

@ -0,0 +1,70 @@
//==============================================================================
//!
//! \file LagrangeFields3D.h
//!
//! \date Jun 8 2011
//!
//! \author Runar Holdahl / SINTEF
//!
//! \brief Class for Lagrange-based finite element vector fields in 3D.
//!
//==============================================================================
#ifndef _LAGRANGE_FIELDS_3D_H
#define _LAGRANGE_FIELDS_3D_H
#include "Fields.h"
/*!
\brief Class for Lagrange-based finite element vector fields in 3D.
\details This class implements the functions required to evaluate a 3D
Lagrange vector field at a given point in parametrical or physical coordinates.
*/
class LagrangeFields3D : public Fields
{
public:
//! \brief The constructor sets the field name.
//! \param[in] geometry Spline volume geometry
//! \param[in] name Name of spline field
LagrangeFields3D(Matrix X, int nx, int ny, int nz, int px, int py, int pz, char* name = NULL);
//! \brief Empty destructor.
virtual ~LagrangeFields3D() {}
// Methods to compute field values
//================================
//! \brief Computes the value in a given node/control point.
//! \param[in] node Node number
//! \param[out] vals Node values
bool valueNode(int node, Vector& vals) const;
//! \brief Computes the value at a given local coordinate.
//! \param[in] fe Finite element definition
//! \param[out] vals Values in local point in given element
bool valueFE(const FiniteElement& fe, Vector& vals) const;
//! \brief Computed the value at a given global coordinate.
//! \param[in] x Global/physical coordinate for point
//! \param[in] vals Values in given physical coordinate
bool valueCoor(const Vec3& x, Vector& vals) const;
//! \brief Computes the gradient for a given local coordinate.
//! \param[in] fe Finite element
//! \param[out] grad Gradient of solution in a given local coordinate
bool gradFE(const FiniteElement& fe, Matrix& grad) const;
//! \brief Computes the gradient for a given global/physical coordinate.
//! \param[in] x Global coordinate
//! \param[out] grad Gradient of solution in a given global coordinate
bool gradCoor(const Vec3& x, Matrix& grad) const;
protected:
Matrix coord; //!< Matrix of nodel coordinates
int n1, n2, n3; //!< Number of nodes in each parameter direction
int p1, p2, p3; //!< Element order in each parameter direction
};
#endif