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

Added: Calculation of characteristic element size in the FE framework

Browse Source
This commit is contained in:
Knut Morten Okstad 2017-04-26 13:53:38 +02:00
parent 21eb5a7bc3
commit 0b0a365b38
17 changed files with 127 additions and 101 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
src/ASM/ASM2D.C
View File

@ -22,6 +22,7 @@
#include "LR/ASMu2DIB.h"
#include "LR/ASMu2Dmx.h"
#endif
#include "Vec3Oper.h"
ASMbase* ASM2D::create (ASM::Discretization discretization, unsigned char nf)
@ -101,3 +102,10 @@ ASMbase* ASM2D::clone (const CharVec& nf) const
#undef TRY_CLONE1
#undef TRY_CLONE2
double ASM2D::getElementSize (const std::vector<Vec3>& XC)
{
// Find longest diagonal (diameter av minste omskrivende sirkel)
return std::max((XC[3]-XC[0]).length(),(XC[2]-XC[1]).length());
}

9
src/ASM/ASM2D.h
View File

@ -19,6 +19,7 @@
#include <cstddef>
class ASMbase;
class Vec3;
/*!
@ -63,7 +64,7 @@ public:
ASMbase* clone(const CharVec& nf = CharVec()) const;
//! \brief Checks that the patch is modelled in a right-hand-side system.
virtual bool checkRightHandSystem() = 0;
virtual bool checkRightHandSystem() { return false; }
//! \brief Refines the parametrization by inserting extra knots uniformly.
//! \param[in] dir Parameter direction to refine
@ -129,7 +130,11 @@ public:
int dir, int nSegSpan) const = 0;
//! \brief Returns the node index for a given corner.
virtual int getCorner(int I, int J, int basis = 1) const = 0;
virtual int getCorner(int I, int J, int basis) const = 0;
protected:
//! \brief Returns characteristic element size based on corner coordinates.
static double getElementSize(const std::vector<Vec3>& XC);
};
#endif

13
src/ASM/ASM3D.C
View File

@ -14,11 +14,11 @@
#include "ASM3D.h"
#include "ASMs3DSpec.h"
#include "ASMs3Dmx.h"
#include "ASMs3DLag.h"
#include "ASMs3DmxLag.h"
#ifdef HAS_LRSPLINE
#include "LR/ASMu3D.h"
#endif
#include "Vec3Oper.h"
ASMbase* ASM3D::create (ASM::Discretization discretization, unsigned char nf)
@ -81,3 +81,14 @@ ASMbase* ASM3D::clone (const CharVec& nf) const
#undef TRY_CLONE1
#undef TRY_CLONE2
double ASM3D::getElementSize (const std::vector<Vec3>& XC)
{
// Find the longest diagonal (diameter av minste omvskrivende kule)
double siz = (XC[7] - XC[0]).length();
for (int c = 1; c < 4; c++)
siz = std::max(siz,(XC[7-c]-XC[c]).length());
return siz;
}

12
src/ASM/ASM3D.h
View File

@ -19,6 +19,7 @@
#include <cstddef>
class ASMbase;
class Vec3;
/*!
@ -63,7 +64,7 @@ public:
ASMbase* clone(const CharVec& nf = CharVec()) const;
//! \brief Checks that the patch is modelled in a right-hand-side system.
virtual bool checkRightHandSystem() = 0;
virtual bool checkRightHandSystem() { return false; }
//! \brief Refines the parametrization by inserting extra knots uniformly.
//! \param[in] dir Parameter direction to refine
@ -161,10 +162,13 @@ public:
int dir, int nSegSpan) const = 0;
//! \brief Returns the node index for a given corner.
virtual int getCorner(int I, int J, int K, int basis = 1) const = 0;
virtual int getCorner(int I, int J, int K, int basis) const = 0;
//! \brief Returns the node indices for a given edge.
virtual std::vector<int> getEdge(int lEdge, bool open, int basis = 1) const = 0;
virtual std::vector<int> getEdge(int lEdge, bool open, int basis) const = 0;
protected:
//! \brief Returns characteristic element size based on corner coordinates.
static double getElementSize(const std::vector<Vec3>& XC);
};
#endif

11
src/ASM/ASMs1D.C
View File

@ -667,7 +667,7 @@ const Vector& ASMs1D::getGaussPointParameters (Matrix& uGP, int nGauss,
}
void ASMs1D::getElementEnds (int i, Vec3Vec& XC) const
double ASMs1D::getElementEnds (int i, Vec3Vec& XC) const
{
RealArray::const_iterator uit = curv->basis().begin();
@ -687,11 +687,14 @@ void ASMs1D::getElementEnds (int i, Vec3Vec& XC) const
for (int j = 0; j < 2; j++, pt += dim)
XC.push_back(Vec3(pt,dim));
if (elmCS.empty()) return;
// Calculate the element length
double h = (XC.back() - XC.front()).length();
if (elmCS.empty()) return h;
// Add the local Z-axis as the third vector
int iel = i - curv->order();
XC.push_back(elmCS[iel][2]);
return h;
}
@ -800,7 +803,7 @@ bool ASMs1D::integrate (Integrand& integrand,
if (!this->getElementCoordinates(fe.Xn,1+iel)) return false;
if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS)
this->getElementEnds(p1+iel,fe.XC);
fe.h = this->getElementEnds(p1+iel,fe.XC);
if (integrand.getIntegrandType() & Integrand::NODAL_ROTATIONS)
{
@ -951,7 +954,7 @@ bool ASMs1D::integrate (Integrand& integrand, int lIndex,
if (!this->getElementCoordinates(fe.Xn,1+iel)) return false;
if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS)
this->getElementEnds(iel+curv->order(),fe.XC);
fe.h = this->getElementEnds(iel+curv->order(),fe.XC);
if (integrand.getIntegrandType() & Integrand::NODAL_ROTATIONS)
{

3
src/ASM/ASMs1D.h
View File

@ -325,7 +325,8 @@ protected:
//! \brief Computes the element end coordinates.
//! \param[in] i Parameter index for the knot-span element
//! \param[out] XC Coordinates of the element corners
void getElementEnds(int i, std::vector<Vec3>& XC) const;
//! \return Element length
double getElementEnds(int i, std::vector<Vec3>& XC) const;
//! \brief Returns nodal rotation matrices for an element, if any.
//! \param[out] T Array of nodal rotation matrices

8
src/ASM/ASMs1DLag.C
View File

@ -179,7 +179,7 @@ bool ASMs1DLag::updateCoords (const Vector& displ)
\brief Extracts the element end points from the element coordinates.
*/
static void getEndPoints (const Matrix& Xnod, Vec3Vec& XC)
static double getEndPoints (const Matrix& Xnod, Vec3Vec& XC)
{
XC.resize(2);
for (size_t i = 0; i < Xnod.cols(); i++)
@ -187,6 +187,8 @@ static void getEndPoints (const Matrix& Xnod, Vec3Vec& XC)
XC[0][i] = Xnod(i+1,1);
XC[1][i] = Xnod(i+1,Xnod.rows());
}
return (XC.back() - XC.front()).length(); // Element length
}
@ -236,7 +238,7 @@ bool ASMs1DLag::integrate (Integrand& integrand,
if (!this->getElementCoordinates(fe.Xn,1+iel)) return false;
if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS)
getEndPoints(fe.Xn,fe.XC);
fe.h = getEndPoints(fe.Xn,fe.XC);
if (integrand.getIntegrandType() & Integrand::NODAL_ROTATIONS)
{
@ -358,7 +360,7 @@ bool ASMs1DLag::integrate (Integrand& integrand, int lIndex,
if (!this->getElementCoordinates(fe.Xn,1+iel)) return false;
if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS)
getEndPoints(fe.Xn,fe.XC);
fe.h = getEndPoints(fe.Xn,fe.XC);
if (integrand.getIntegrandType() & Integrand::NODAL_ROTATIONS)
{

14
src/ASM/ASMs2D.C
View File

@ -1384,7 +1384,7 @@ void ASMs2D::getElementBorders (int i1, int i2, double* u, double* v) const
}
void ASMs2D::getElementCorners (int i1, int i2, Vec3Vec& XC) const
double ASMs2D::getElementCorners (int i1, int i2, Vec3Vec& XC) const
{
// Fetch parameter values of the element (knot-span) corners
RealArray u(2), v(2);
@ -1401,6 +1401,8 @@ void ASMs2D::getElementCorners (int i1, int i2, Vec3Vec& XC) const
const double* pt = &XYZ.front();
for (int i = 0; i < 4; i++, pt += dim)
XC.push_back(Vec3(pt,dim));
return getElementSize(XC);
}
@ -1520,7 +1522,7 @@ bool ASMs2D::integrate (Integrand& integrand,
}
if (useElmVtx)
this->getElementCorners(i1-1,i2-1,fe.XC);
fe.h = this->getElementCorners(i1-1,i2-1,fe.XC);
if (integrand.getIntegrandType() & Integrand::G_MATRIX)
{
@ -1787,11 +1789,11 @@ bool ASMs2D::integrate (Integrand& integrand,
if (integrand.getIntegrandType() & Integrand::ELEMENT_CENTER)
{
// Compute the element center
this->getElementCorners(i1-1,i2-1,fe.XC);
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]);
}
else if (useElmVtx)
this->getElementCorners(i1-1,i2-1,fe.XC);
fe.h = this->getElementCorners(i1-1,i2-1,fe.XC);
if (integrand.getIntegrandType() & Integrand::G_MATRIX)
{
@ -1936,7 +1938,7 @@ bool ASMs2D::integrate (Integrand& integrand,
this->getElementBorders(i1-1,i2-1,u,v);
if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS)
this->getElementCorners(i1-1,i2-1,fe.XC);
fe.h = this->getElementCorners(i1-1,i2-1,fe.XC);
// Initialize element quantities
LocalIntegral* A = integrand.getLocalIntegral(MNPC[jel].size(),fe.iel);
@ -2141,7 +2143,7 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex,
if (!this->getElementCoordinates(Xnod,iel)) return false;
if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS)
this->getElementCorners(i1-1,i2-1,fe.XC);
fe.h = this->getElementCorners(i1-1,i2-1,fe.XC);
if (integrand.getIntegrandType() & Integrand::G_MATRIX)
{

5
src/ASM/ASMs2D.h
View File

@ -190,7 +190,7 @@ public:
bool local = false) const;
//! \brief Returns the node index for a given corner.
virtual int getCorner(int I, int J, int basis = 1) const;
virtual int getCorner(int I, int J, int basis) const;
//! \brief Assigns new global node numbers for all nodes of the patch.
//! \param nodes Object with global nodes numbers to assign to this patch
@ -527,7 +527,8 @@ protected:
//! \param[in] i1 Parameter index in u-direction
//! \param[in] i2 Parameter index in v-direction
//! \param[out] XC Coordinates of the element corners
void getElementCorners(int i1, int i2, std::vector<Vec3>& XC) const;
//! \return Characteristic element size
double getElementCorners(int i1, int i2, std::vector<Vec3>& XC) const;
using ASMbase::generateThreadGroups;
//! \brief Generates element groups for multi-threading of interior integrals.

12
src/ASM/ASMs3D.C
View File

@ -1640,7 +1640,7 @@ void ASMs3D::getElementBorders (int i1, int i2, int i3,
}
void ASMs3D::getElementCorners (int i1, int i2, int i3, Vec3Vec& XC) const
double ASMs3D::getElementCorners (int i1, int i2, int i3, Vec3Vec& XC) const
{
// Fetch parameter values of the element (knot-span) corners
RealArray u(2), v(2), w(2);
@ -1657,6 +1657,8 @@ void ASMs3D::getElementCorners (int i1, int i2, int i3, Vec3Vec& XC) const
const double* pt = &XYZ.front();
for (int i = 0; i < 8; i++, pt += dim)
XC.push_back(Vec3(pt,dim));
return getElementSize(XC);
}
@ -1767,7 +1769,7 @@ bool ASMs3D::integrate (Integrand& integrand,
}
if (useElmVtx)
this->getElementCorners(i1-1,i2-1,i3-1,fe.XC);
fe.h = this->getElementCorners(i1-1,i2-1,i3-1,fe.XC);
if (integrand.getIntegrandType() & Integrand::G_MATRIX)
{
@ -2048,12 +2050,12 @@ bool ASMs3D::integrate (Integrand& integrand,
if (integrand.getIntegrandType() & Integrand::ELEMENT_CENTER)
{
// Compute the element center
this->getElementCorners(i1-1,i2-1,i3-1,fe.XC);
fe.h = this->getElementCorners(i1-1,i2-1,i3-1,fe.XC);
X = 0.125*(fe.XC[0]+fe.XC[1]+fe.XC[2]+fe.XC[3]+
fe.XC[4]+fe.XC[5]+fe.XC[6]+fe.XC[7]);
}
else if (useElmVtx)
this->getElementCorners(i1-1,i2-1,i3-1,fe.XC);
fe.h = this->getElementCorners(i1-1,i2-1,i3-1,fe.XC);
if (integrand.getIntegrandType() & Integrand::G_MATRIX)
{
@ -2284,7 +2286,7 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex,
}
if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS)
this->getElementCorners(i1-1,i2-1,i3-1,fe.XC);
fe.h = this->getElementCorners(i1-1,i2-1,i3-1,fe.XC);
if (integrand.getIntegrandType() & Integrand::G_MATRIX)
{

8
src/ASM/ASMs3D.h
View File

@ -206,10 +206,9 @@ public:
bool local = false) const;
//! \brief Returns the node index for a given corner.
virtual int getCorner(int I, int J, int K, int basis = 1) const;
virtual int getCorner(int I, int J, int K, int basis) const;
//! \brief Returns the node indices for a given edge.
virtual std::vector<int> getEdge(int lEdge, bool open, int basis = 1) const;
virtual std::vector<int> getEdge(int lEdge, bool open, int basis) const;
//! \brief Assigns new global node numbers for all nodes of the patch.
//! \param nodes Object with global nodes numbers to assign to this patch
@ -595,7 +594,8 @@ protected:
//! \param[in] i2 Parameter index in v-direction
//! \param[in] i3 Parameter index in w-direction
//! \param[out] XC Coordinates of the element corners
void getElementCorners(int i1, int i2, int i3, std::vector<Vec3>& XC) const;
//! \return Characteristic element size
double getElementCorners(int i1, int i2, int i3, std::vector<Vec3>& XC) const;
//! \brief Generates element groups for multi-threading of interior integrals.
//! \param[in] integrand Object with problem-specific data and methods

8
src/ASM/FiniteElement.C
View File

@ -14,18 +14,12 @@
#include "FiniteElement.h"
std::ostream& operator<< (std::ostream& os, const FiniteElement& fe)
{
return fe.write(os);
}
std::ostream& FiniteElement::write (std::ostream& os) const
{
os <<"FiniteElement: iel="<< iel <<" iGP="<< iGP <<" p="<< p
<<"\n u, v, w: "<< u <<" "<< v <<" "<< w
<<"\n xi, eta, zeta: "<< xi <<" "<< eta <<" "<< zeta
<<"\n detJxW: "<< detJxW << std::endl;
<<"\n h, detJxW: "<< h <<" "<< detJxW << std::endl;
if (!N.empty()) os <<"N:"<< N;
if (!dNdX.empty()) os <<"dNdX:"<< dNdX;
if (!d2NdX2.empty()) os <<"d2NdX2:"<< d2NdX2;

22
src/ASM/FiniteElement.h
View File

@ -27,8 +27,8 @@ class FiniteElement
{
public:
//! \brief Default constructor.
explicit FiniteElement(size_t n = 0, size_t i = 0) : iGP(i), N(n), iel(0), p(0), Te(3)
{ u = v = w = xi = eta = zeta = 0.0; detJxW = 1.0; }
explicit FiniteElement(size_t n = 0, size_t i = 0) : iGP(i), N(n), Te(3)
{ iel = p = 0; u = v = w = xi = eta = zeta = h = 0.0; detJxW = 1.0; }
//! \brief Empty destructor.
virtual ~FiniteElement() {}
@ -52,8 +52,11 @@ protected:
//! \brief Writes the finite element object to the given output stream.
virtual std::ostream& write(std::ostream& os) const;
//! \brief Global Output stream operator.
friend std::ostream& operator<<(std::ostream& os, const FiniteElement& fe);
//! \brief Global output stream operator.
friend std::ostream& operator<<(std::ostream& os, const FiniteElement& fe)
{
return fe.write(os);
}
public:
// Gauss point quantities
@ -68,11 +71,12 @@ public:
Vector N; //!< Basis function values
Matrix dNdX; //!< First derivatives (gradient) of the basis functions
Matrix3D d2NdX2; //!< Second derivatives of the basis functions
Matrix G; //!< Matrix used for stabilized methods
Matrix G; //!< Matrix used for stabilized methods
// Element quantities
int iel; //!< Element identifier
short int p; //!< Polynomial order of the basis functions
double h; //!< Characteristic element size/diameter
Vec3Vec XC; //!< Array with element corner coordinate vectors
Vector Navg; //!< Volume-averaged basis function values
Matrix Xn; //!< Matrix of element nodal coordinates
@ -97,7 +101,7 @@ public:
//! \brief Returns a const reference to the basis function values.
virtual const Vector& basis(char b) const { return b == 1 ? N : Nx[b-2]; }
//! \brief Returns a reference to the basis function values.
Vector& basis(char b) { return b == 1 ? N : Nx[b-2]; }
virtual Vector& basis(char b) { return b == 1 ? N : Nx[b-2]; }
//! \brief Returns a const reference to the basis function derivatives.
virtual const Matrix& grad(char b) const { return b == 1 ? dNdX : dNxdX[b-2]; }
@ -114,9 +118,9 @@ protected:
virtual std::ostream& write(std::ostream& os) const;
private:
Vectors Nx;
std::vector<Matrix> dNxdX;
std::vector<Matrix3D> d2NxdX2;
std::vector<Vector> Nx; //!< Basis function values
std::vector<Matrix> dNxdX; //!< First derivatives of the basis functions
std::vector<Matrix3D> d2NxdX2; //!< Second derivatives of the basis functions
};
#endif

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

@ -60,7 +60,8 @@ bool ASMu2D::read (std::istream& is)
{
if (shareFE) return false;
// read inputfile as either an LRSpline file directly or a tensor product B-spline and convert
// read inputfile as either an LRSpline file directly
// or a tensor product B-spline and convert
char firstline[256];
is.getline(firstline, 256);
if (strncmp(firstline, "# LRSPLINE", 10) == 0) {
@ -135,13 +136,6 @@ void ASMu2D::clear (bool retainGeometry)
}
bool ASMu2D::checkRightHandSystem ()
{
std::cout <<"ASMu2D::checkRightHandSystem(): Not available for LR-splines (ignored)."<< std::endl;
return false;
}
bool ASMu2D::cornerRefine (int minBasisfunctions)
{
if (!lrspline) return false;
@ -930,14 +924,14 @@ void ASMu2D::getGaussPointParameters (RealArray& uGP, int dir, int nGauss,
}
void ASMu2D::getElementCorners (int iel, Vec3Vec& XC) const
double ASMu2D::getElementCorners (int iel, Vec3Vec& XC) const
{
#ifdef INDEX_CHECK
if (iel < 1 || iel > lrspline->nElements())
{
std::cerr <<" *** ASMu2D::getElementCorners: Element index "<< iel
<<" out of range [1,"<< lrspline->nElements() <<"]."<< std::endl;
return;
<<" out of range [1,"<< lrspline->nElements() <<"]."<< std::endl;
return 0.0;
}
#endif
@ -954,6 +948,8 @@ void ASMu2D::getElementCorners (int iel, Vec3Vec& XC) const
lrspline->point(point,u[i],v[i],iel-1);
XC.push_back(SplineUtils::toVec3(point,nsd));
}
return getElementSize(XC);
}
@ -1031,7 +1027,7 @@ bool ASMu2D::integrate (Integrand& integrand,
}
if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS)
this->getElementCorners(iel,fe.XC);
fe.h = this->getElementCorners(iel,fe.XC);
if (integrand.getIntegrandType() & Integrand::ELEMENT_CENTER)
{
@ -1054,7 +1050,7 @@ bool ASMu2D::integrate (Integrand& integrand,
if (xr)
{
// --- Selective reduced integration loop --------------------------------
// --- Selective reduced integration loop ------------------------------
for (int j = 0; j < nRed; j++)
for (int i = 0; i < nRed; i++)
@ -1089,7 +1085,7 @@ bool ASMu2D::integrate (Integrand& integrand,
}
}
// --- Integration loop over all Gauss points in each direction ------------
// --- Integration loop over all Gauss points in each direction ----------
int jp = (iel-1)*nGauss*nGauss;
fe.iGP = firstIp + jp; // Global integration point counter
@ -1248,12 +1244,12 @@ bool ASMu2D::integrate (Integrand& integrand,
}
if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS)
this->getElementCorners(iel,fe.XC);
fe.h = this->getElementCorners(iel,fe.XC);
if (integrand.getIntegrandType() & Integrand::ELEMENT_CENTER)
{
// Compute the element center
this->getElementCorners(iel,fe.XC);
fe.h = this->getElementCorners(iel,fe.XC);
X = 0.25*(fe.XC[0]+fe.XC[1]+fe.XC[2]+fe.XC[3]);
}
@ -1265,7 +1261,7 @@ bool ASMu2D::integrate (Integrand& integrand,
continue;
}
// --- Integration loop over all quadrature points in this element ---------
// --- Integration loop over all quadrature points in this element -------
size_t jp = MPitg[iel-1]; // Patch-wise integration point counter
fe.iGP = firstIp + jp; // Global integration point counter
@ -1423,9 +1419,9 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex,
this->getGaussPointParameters(gpar[t2-1],t2-1,nGP,iel,xg);
if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS)
this->getElementCorners(iel,fe.XC);
fe.h = this->getElementCorners(iel,fe.XC);
// --- Integration loop over all Gauss points along the edge -------------
// --- Integration loop over all Gauss points along the edge ---------------
fe.iGP = firstp; // Global integration point counter
firstp += nGP;

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

@ -109,9 +109,6 @@ public:
//! \brief Returns the total number of nodes in this patch.
virtual size_t getNoNodes(int basis = 0) const;
//! \brief Checks that the patch is modelled in a right-hand-side system.
virtual bool checkRightHandSystem();
//! \brief Refines along the diagonal of the LR-spline patch.
//! \details Progressively refine until the LR-spline object contains at least
//! \a minBasisfunctions basis functions.
@ -160,7 +157,8 @@ public:
//! \param[in] open If \e true, exclude the end points of the edge
//! \param[in] dof Which DOFs to constrain at each node on the edge
//! \param[in] code Inhomogeneous dirichlet condition code
virtual void constrainEdge(int dir, bool open, int dof, int code, char);
//! \param[in] basis Which basis to constrain edge for
virtual void constrainEdge(int dir, bool open, int dof, int code, char basis);
//! \brief Constrains all DOFs in local directions on a given boundary edge.
//! \param[in] dir Parameter direction defining the edge to constrain
//! \param[in] open If \e true, exclude the end points of the edge
@ -242,9 +240,10 @@ public:
const std::map<int,int>* g2l = nullptr);
//! \brief Returns the node index for a given corner.
virtual int getCorner(int I, int J, int basis = 1) const;
virtual int getCorner(int I, int J, int basis) const;
//! \brief Returns the node indices for a given edge.
std::vector<int> getEdgeNodes(int edge, int basis = 1) const;
std::vector<int> getEdgeNodes(int edge, int basis) const;
protected:
//! \brief Evaluates an integral over the interior patch domain.
@ -468,7 +467,8 @@ protected:
//! \brief Computes the element corner coordinates.
//! \param[in] iel Element index
//! \param[out] XC Coordinates of the element corners
void getElementCorners(int iel, std::vector<Vec3>& XC) const;
//! \return Characteristic element size
double getElementCorners(int iel, std::vector<Vec3>& XC) const;
//! \brief Evaluates the basis functions and derivatives of order \a derivs
//! of an element.
@ -503,8 +503,8 @@ protected:
const std::vector<Matrix>& bezierExtract; //!< Bezier extraction matrices
std::vector<Matrix> myBezierExtract; //!< Bezier extraction matrices
Go::BsplineBasis bezier_u;
Go::BsplineBasis bezier_v;
Go::BsplineBasis bezier_u; //!< Bezier basis in the u-direction
Go::BsplineBasis bezier_v; //!< Bezier basis in the v-direction
ThreadGroups threadGroups; //!< Element groups for multi-threaded assembly
};

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

@ -124,13 +124,6 @@ void ASMu3D::clear (bool retainGeometry)
}
bool ASMu3D::checkRightHandSystem ()
{
std::cout <<"ASMu3D::checkRightHandSystem(): Not available for LR-splines (ignored)."<< std::endl;
return false;
}
bool ASMu3D::refine (int dir, const RealArray& xi)
{
if (!tensorspline || dir < 0 || dir > 2 || xi.empty()) return false;
@ -795,7 +788,7 @@ void ASMu3D::getGaussPointParameters (RealArray& uGP, int dir, int nGauss,
}
void ASMu3D::getElementCorners (int iEl, Vec3Vec& XC) const
double ASMu3D::getElementCorners (int iEl, Vec3Vec& XC) const
{
LR::Element* el = lrspline->getElement(iEl);
double u[2] = { el->getParmin(0), el->getParmax(0) };
@ -813,6 +806,8 @@ void ASMu3D::getElementCorners (int iEl, Vec3Vec& XC) const
lrspline->point(pt,u[i],v[j],w[k],iEl);
XC.push_back(SplineUtils::toVec3(pt));
}
return getElementSize(XC);
}
@ -1020,7 +1015,7 @@ bool ASMu3D::integrate (Integrand& integrand,
if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS)
this->getElementCorners(iEl, fe.XC);
fe.h = this->getElementCorners(iEl,fe.XC);
if (integrand.getIntegrandType() & Integrand::G_MATRIX)
{
@ -1031,7 +1026,7 @@ bool ASMu3D::integrate (Integrand& integrand,
}
else if (integrand.getIntegrandType() & Integrand::AVERAGE)
{
// --- Compute average value of basis functions over the element -----
// --- Compute average value of basis functions over the element -------
fe.Navg.resize(nBasis,true);
double vol = 0.0;
@ -1081,7 +1076,7 @@ bool ASMu3D::integrate (Integrand& integrand,
std::cerr << "Haven't really figured out what this part does yet\n";
exit(42142);
#if 0
// --- Selective reduced integration loop ----------------------------
// --- Selective reduced integration loop ------------------------------
int ip = (((i3-p3)*nRed*nel2 + i2-p2)*nRed*nel1 + i1-p1)*nRed;
for (int k = 0; k < nRed; k++, ip += nRed*(nel2-1)*nRed*nel1)
@ -1117,7 +1112,7 @@ bool ASMu3D::integrate (Integrand& integrand,
}
// --- Integration loop over all Gauss points in each direction --------
// --- Integration loop over all Gauss points in each direction ----------
fe.iGP = iEl*nGauss*nGauss*nGauss; // Global integration point counter
@ -1331,7 +1326,7 @@ bool ASMu3D::integrate (Integrand& integrand, int lIndex,
}
if (integrand.getIntegrandType() & Integrand::ELEMENT_CORNERS)
this->getElementCorners(iEl,fe.XC);
fe.h = this->getElementCorners(iEl,fe.XC);
if (integrand.getIntegrandType() & Integrand::G_MATRIX)
{
@ -1350,7 +1345,7 @@ bool ASMu3D::integrate (Integrand& integrand, int lIndex,
break;
}
// --- Integration loop over all Gauss points in each direction --------
// --- Integration loop over all Gauss points in each direction ------------
fe.iGP = firstp; // Global integration point counter
int k1,k2,k3;
@ -1560,7 +1555,7 @@ bool ASMu3D::integrateEdge (Integrand& integrand, int lEdge,
bool ok = integrand.initElementBou(MNPC[iel-1],*A);
// --- Integration loop over all Gauss points along the edge -----------
// --- Integration loop over all Gauss points along the edge -----------------
fe.iGP = firstp + ip; // Global integration point counter

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

@ -89,17 +89,17 @@ public:
std::vector<int> getFaceNodes(int face, int basis = 1) const;
//! \brief Finds the global (or patch-local) node numbers on a patch boundary.
//! \param[in] lIndex Local index of the boundary edge
//! \param glbNodes Array of global boundary node numbers
//! \param local If \e true return patch-local node numbers
//! \param[in] lIndex Local index of the boundary face
//! \param nodes Array of node numbers
//! \param basis Which basis to grab nodes for (0 for all)
//! \param local If \e true, return patch-local node numbers
virtual void getBoundaryNodes(int lIndex, IntVec& nodes,
int basis, int, bool local) const;
//! \brief Returns the node index for a given corner.
virtual int getCorner(int I, int J, int K, int basis = 1) const;
virtual int getCorner(int I, int J, int K, int basis) const;
//! \brief Returns the node indices for a given edge.
virtual std::vector<int> getEdge(int lEdge, bool open, int basis = 1) const;
virtual std::vector<int> getEdge(int lEdge, bool open, int basis) const;
//! \brief Returns the polynomial order in each parameter direction.
//! \param[out] p1 Order in first (u) direction
@ -107,9 +107,6 @@ public:
//! \param[out] p3 Order in third (w) direction
virtual bool getOrder(int& p1, int& p2, int& p3) const;
//! \brief Checks that the patch is modelled in a right-hand-side system.
virtual bool checkRightHandSystem();
//! \brief Refines the parametrization by inserting tensor knots uniformly.
//! \param[in] dir Parameter direction to refine
//! \param[in] nInsert Number of extra knots to insert in each knot-span
@ -418,7 +415,8 @@ protected:
//! \brief Computes the element corner coordinates.
//! \param[in] iel Element index
//! \param[out] XC Coordinates of the element corners
void getElementCorners(int iel, std::vector<Vec3>& XC) const;
//! \return Characteristic element size
double getElementCorners(int iel, std::vector<Vec3>& XC) const;
//! \brief Evaluate all basis functions and \a derivs number of derivatives on one element
virtual void evaluateBasis(FiniteElement &el, int derivs) const;