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

Added: Driver for integration of global nodal forces on boundaries

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@2359 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo
2013-05-07 17:01:48 +00:00
committed by Knut Morten Okstad
parent 58bf05a58b
commit f51db0a2ac
13 changed files with 593 additions and 768 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
src/ASM/ASMbase.h
View File

@@ -174,6 +174,11 @@ public:
//! in one element
virtual bool getElementCoordinates(Matrix& X, int iel) const = 0;
//! \brief Finds the global numbers of the nodes on a patch boundary.
//! \param[in] lIndex Local index of the boundary face/edge
//! \param glbNodes Array of global boundary node numbers
virtual void getBoundaryNodes(int lIndex, IntVec& glbNodes) const = 0;
//! \brief Prints out the nodal coordinates of this patch to the given stream.
void printNodes(std::ostream& os, const char* heading = 0) const;

15
src/ASM/ASMs1D.C
View File

@@ -452,6 +452,21 @@ bool ASMs1D::updateCoords (const Vector& displ)
}
void ASMs1D::getBoundaryNodes (int lIndex, IntVec& glbNodes) const
{
if (!curv) return; // silently ignore empty patches
int iel = lIndex == 1 ? 0 : this->getNoElms()-1;
if (MLGE[iel] > 0)
{
if (lIndex == 1)
glbNodes.push_back(MLGN[MNPC[iel].back()]);
else if (lIndex == 2)
glbNodes.push_back(MLGN[MNPC[iel].front()]);
}
}
bool ASMs1D::getOrder (int& p1, int& p2, int& p3) const
{
p2 = p3 = 0;

8
src/ASM/ASMs1D.h
View File

@@ -36,9 +36,10 @@ public:
//! \brief Empty destructor.
virtual ~ASMs1D() {}
//! \brief Returns the spline surface representing the geometry of this patch.
//! \brief Returns the spline curve representing the geometry of this patch.
Go::SplineCurve* getCurve() const { return curv; }
// Methods for model generation
// ============================
@@ -76,6 +77,11 @@ public:
//! \param[in] displ Incremental displacements to update the coordinates with
virtual bool updateCoords(const Vector& displ);
//! \brief Finds the global number of the node on a patch end.
//! \param[in] lIndex Local index of the end point
//! \param glbNodes Array of global boundary node numbers
virtual void getBoundaryNodes(int lIndex, IntVec& glbNodes) const;
//! \brief Refines the parametrization by inserting extra knots.
//! \param[in] xi Relative positions of added knots in each existing knot span
bool refine(const RealArray& xi);

50
src/ASM/ASMs2D.C
View File

@@ -1181,6 +1181,56 @@ bool ASMs2D::updateCoords (const Vector& displ)
}
void ASMs2D::getBoundaryNodes (int lIndex, IntVec& nodes) const
{
if (!surf) return; // silently ignore empty patches
const int p1 = surf->order_u();
const int p2 = surf->order_v();
const int n1 = surf->numCoefs_u();
const int n2 = surf->numCoefs_v();
#if SP_DEBUG > 1
size_t last = nodes.size();
#endif
int iel = 0;
for (int i2 = p2; i2 <= n2; i2++)
for (int i1 = p1; i1 <= n1; i1++, iel++)
if (surf->knotSpan(0,i1-1) > 0.0)
if (surf->knotSpan(1,i2-1) > 0.0)
{
int inod = 0, lnod = 0;
for (int j2 = p2; j2 > 0; j2--)
for (int j1 = p1; j1 > 0; j1--, lnod++)
{
if (lIndex == 1 && i1 == p1 && j1 == p1)
inod = MNPC[iel][lnod]; // left edge
else if (lIndex == 3 && i2 == p2 && j2 == p2)
inod = MNPC[iel][lnod]; // bottom edge
else if (lIndex == 2 && i1 == n1 && j1 == 1)
inod = MNPC[iel][lnod]; // right edge
else if (lIndex == 4 && i2 == n2 && j2 == 1)
inod = MNPC[iel][lnod]; // top edge
else
continue;
inod = MLGN[inod];
if (std::find(nodes.begin(),nodes.end(),inod) == nodes.end())
nodes.push_back(inod);
}
}
#if SP_DEBUG > 1
std::cout <<"Boundary nodes in patch "<< idx+1 <<" edge "<< lIndex <<":";
if (nodes.size() == last)
std::cout <<" (none)";
else for (size_t i = last; i < nodes.size(); i++)
std::cout <<" "<< nodes[i];
std::cout << std::endl;
#endif
}
bool ASMs2D::getOrder (int& p1, int& p2) const
{
if (!surf) return false;

5
src/ASM/ASMs2D.h
View File

@@ -161,6 +161,11 @@ public:
//! \param[in] displ Incremental displacements to update the coordinates with
virtual bool updateCoords(const Vector& displ);
//! \brief Finds the global numbers of the nodes on a patch boundary.
//! \param[in] lIndex Local index of the boundary edge
//! \param glbNodes Array of global boundary node numbers
virtual void getBoundaryNodes(int lIndex, IntVec& glbNodes) 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
//! \param[in] basis Which basis to assign node numbers for in mixed methods

6
src/ASM/ASMs3D.C
View File

@@ -1487,6 +1487,12 @@ bool ASMs3D::updateCoords (const Vector& displ)
}
void ASMs3D::getBoundaryNodes (int lIndex, IntVec& nodes) const
{
// TODO: Implement this before attempting 3D FSI simulations
}
bool ASMs3D::getOrder (int& p1, int& p2, int& p3) const
{
if (!svol) return false;

5
src/ASM/ASMs3D.h
View File

@@ -180,6 +180,11 @@ public:
//! \param[in] displ Incremental displacements to update the coordinates with
virtual bool updateCoords(const Vector& displ);
//! \brief Finds the global numbers of the nodes on a patch boundary.
//! \param[in] lIndex Local index of the boundary face
//! \param glbNodes Array of global boundary node numbers
virtual void getBoundaryNodes(int lIndex, IntVec& glbNodes) 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
//! \param[in] basis Which basis to assign node numbers for in mixed methods

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

@@ -291,7 +291,6 @@ bool ASMu2D::raiseOrder (int ru, int rv)
}
bool ASMu2D::generateFEMTopology ()
{
// At this point we are through with the tensor spline object.
@@ -1487,3 +1486,9 @@ bool ASMu2D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
return true;
}
void ASMu2D::getBoundaryNodes (int lIndex, IntVec& nodes) const
{
// TODO: Implement this before attempting FSI simulations with LR B-splines
}

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

@@ -82,6 +82,11 @@ public:
//! \param[in] displ Incremental displacements to update the coordinates with
virtual bool updateCoords(const Vector& displ);
//! \brief Finds the global numbers of the nodes on a patch boundary.
//! \param[in] lIndex Local index of the boundary edge
//! \param glbNodes Array of global boundary node numbers
virtual void getBoundaryNodes(int lIndex, IntVec& glbNodes) const;
//! \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.

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

@@ -44,57 +44,17 @@
ASMu3D::ASMu3D (unsigned char n_f)
: ASMunstruct(3,3,n_f), lrspline(0), tensorspline(0)
: ASMunstruct(3,3,n_f), lrspline(0), tensorspline(0)
{
ASMunstruct::resetNumbering(); // Replace this when going multi-patch...
ASMunstruct::resetNumbering(); // Replace this when going multi-patch...
}
ASMu3D::ASMu3D (const ASMu3D& patch, unsigned char n_f)
: ASMunstruct(patch,n_f), lrspline(patch.lrspline), tensorspline(0)
: ASMunstruct(patch,n_f), lrspline(patch.lrspline), tensorspline(0)
{
}
size_t ASMu3D::getNodeIndex (int globalNum, bool noAddedNodes) const
{
IntVec::const_iterator it = std::find(MLGN.begin(),MLGN.end(),globalNum);
if (it == MLGN.end()) return 0;
size_t inod = 1 + (it-MLGN.begin());
if (noAddedNodes && !xnMap.empty())
{
std::map<size_t,size_t>::const_iterator it = xnMap.find(inod);
if (it != xnMap.end()) return it->second;
}
return inod;
}
/*
char ASMu3D::getNodeType (size_t inod) const
{
std::cerr << "ASM3D::getNodeType not implemented properly\n";
exit(123213);
if (this->isLMn(inod)) return 'L';
// return inod > nodeInd.size() ? 'X' : 'D';
return 0;
}
*/
LR::LRSplineSurface* ASMu3D::getBoundary (int dir)
{
std::cerr << "ASMu3D::getBoundary not implemented properly yet" << std::endl;
exit(776654);
if (dir < -3 || dir == 0 || dir > 3)
return NULL;
// The boundary surfaces are stored internally in the SplineVolume object
// int iface = dir > 0 ? 2*dir-1 : -2*dir-2;
return NULL;
// return lrspline->getBoundarySurface(iface).get();
}
bool ASMu3D::read (std::istream& is)
{
@@ -166,14 +126,6 @@ void ASMu3D::clear (bool retainGeometry)
// Erase the FE data
this->ASMbase::clear(retainGeometry);
xnMap.clear();
nxMap.clear();
}
size_t ASMu3D::getNoNodes (int basis) const
{
return lrspline->nBasisFunctions();
}
@@ -327,6 +279,9 @@ bool ASMu3D::generateFEMTopology ()
}
/*
// this is connecting multiple patches and handling deformed geometries.
// We'll deal with at a later time, for now we only allow single patch models
bool ASMu3D::connectPatch (int face, ASMu3D& neighbor, int nface, int norient)
{
@@ -337,9 +292,6 @@ bool ASMu3D::connectPatch (int face, ASMu3D& neighbor, int nface, int norient)
bool ASMu3D::connectBasis (int face, ASMu3D& neighbor, int nface, int norient,
int basis, int slave, int master)
{
std::cerr << "ASMu3D::connectBasis(...) is not properly implemented yet :(" << std::endl;
exit(776654);
#if 0
if (shareFE && neighbor.shareFE)
return true;
else if (shareFE || neighbor.shareFE)
@@ -473,16 +425,11 @@ bool ASMu3D::connectBasis (int face, ASMu3D& neighbor, int nface, int norient,
}
return true;
#endif
return false;
}
void ASMu3D::closeFaces (int dir, int basis, int master)
{
std::cerr << "ASMu3D::closeFaces(...) is not properly implemented yet :(" << std::endl;
exit(776654);
#if 0
int n1, n2, n3;
if (basis < 1) basis = 1;
if (!this->getSize(n1,n2,n3,basis)) return;
@@ -507,8 +454,9 @@ void ASMu3D::closeFaces (int dir, int basis, int master)
this->makePeriodic(master,master+n1*n2*(n3-1));
break;
}
#endif
}
*/
/*!
A negative \a code value implies direct evaluation of the Dirichlet condition
@@ -719,78 +667,8 @@ void ASMu3D::constrainNode (double xi, double eta, double zeta,
}
/*!
This method projects the function describing the in-homogeneous Dirichlet
boundary condition onto the spline basis defining the boundary surface,
in order to find the control point values which are used as the prescribed
values of the boundary DOFs.
*/
bool ASMu3D::updateDirichlet (const std::map<int,RealFunc*>& func,
const std::map<int,VecFunc*>& vfunc, double time)
{
// std::cerr << "ASMu3D::updateDirichlet not implemented properly yet\n";
std::cerr << "\nWARNING: ASMu3D::updateDirichlet ignored due to non-projecting boundary condition implementation" << std::endl;
return this->ASMbase::updateDirichlet(func,vfunc,time);
#if 0
std::map<int,RealFunc*>::const_iterator fit;
std::map<int,VecFunc*>::const_iterator vfit;
std::vector<DirichletFace>::const_iterator dit;
std::vector<Ipair>::const_iterator nit;
for (size_t i = 0; i < dirich.size(); i++)
{
// Project the function onto the spline surface basis
Go::SplineSurface* dsurf = 0;
if ((fit = func.find(dirich[i].code)) != func.end())
dsurf = SplineUtils::project(dirich[i].surf,*fit->second,time);
else if ((vfit = vfunc.find(dirich[i].code)) != vfunc.end())
dsurf = SplineUtils::project(dirich[i].surf,*vfit->second,nf,time);
else
{
std::cerr <<" *** ASMu3D::updateDirichlet: Code "<< dirich[i].code
<<" is not associated with any function."<< std::endl;
return false;
}
if (!dsurf)
{
std::cerr <<" *** ASMu3D::updateDirichlet: Projection failure."
<< std::endl;
return false;
}
// Loop over the (interior) nodes (control points) of this boundary surface
for (nit = dirich[i].nodes.begin(); nit != dirich[i].nodes.end(); nit++)
for (int dofs = dirich[i].dof; dofs > 0; dofs /= 10)
{
int dof = dofs%10;
// Find the constraint equation for current (node,dof)
MPC pDOF(MLGN[nit->second-1],dof);
MPCIter mit = mpcs.find(&pDOF);
if (mit == mpcs.end()) continue; // probably a deleted constraint
// Find index to the control point value for this (node,dof) in dsurf
RealArray::const_iterator cit = dsurf->coefs_begin();
if (dsurf->dimension() > 1) // A vector field is specified
cit += (nit->first-1)*dsurf->dimension() + (dof-1);
else // A scalar field is specified at this dof
cit += (nit->first-1);
// Now update the prescribed value in the constraint equation
(*mit)->setSlaveCoeff(*cit);
#if SP_DEBUG > 1
std::cout <<"Updated constraint: "<< **mit;
#endif
}
delete dsurf;
}
// The parent class method takes care of the corner nodes with direct
// evaluation of the Dirichlet functions (since they are interpolatory)
#endif
}
#define DERR -999.99
double ASMu3D::getParametricArea (int iel, int dir) const
{
#ifdef INDEX_CHECK
@@ -822,119 +700,61 @@ double ASMu3D::getParametricArea (int iel, int dir) const
#undef DERR
#if 0
int ASMu3D::coeffInd (size_t inod) const
{
#ifdef INDEX_CHECK
if (inod >= nodeInd.size())
{
std::cerr <<" *** ASMu3D::coeffInd: Node index "<< inod
<<" out of range [0,"<< nodeInd.size() <<">."<< std::endl;
return -1;
}
#endif
const int ni = nodeInd[inod].I;
const int nj = nodeInd[inod].J;
const int nk = nodeInd[inod].K;
return (nk*lrspline->numCoefs(1) + nj)*lrspline->numCoefs(0) + ni;
}
#endif
Vec3 ASMu3D::getCoord (size_t inod) const
{
if (inod <= MLGN.size())
{
// This is a node added due to constraints in local directions.
// Find the corresponding original node (see constrainEdgeLocal)
std::map<size_t,size_t>::const_iterator it = xnMap.find(inod);
if (it != xnMap.end()) inod = it->second;
}
if (inod == 0) return Vec3();
LR::Basisfunction *b = lrspline->getBasisfunction(inod-1);
RealArray::const_iterator cit = b->cp();
return Vec3(*cit,*(cit+1),*(cit+2));
}
bool ASMu3D::getElementCoordinates (Matrix& X, int iel) const
{
#ifdef INDEX_CHECK
if (iel < 1 || (size_t)iel > MNPC.size())
{
std::cerr <<" *** ASMu3D::getElementCoordinates: Element index "<< iel
<<" out of range [1,"<< MNPC.size() <<"]."<< std::endl;
return false;
}
if (iel < 1 || (size_t)iel > MNPC.size())
{
std::cerr <<" *** ASMu3D::getElementCoordinates: Element index "<< iel
<<" out of range [1,"<< MNPC.size() <<"]."<< std::endl;
return false;
}
#endif
const IntVec& mnpc = MNPC[iel-1];
X.resize(3,mnpc.size());
LR::Element* el = lrspline->getElement(iel-1);
X.resize(3,el->nBasisFunctions());
int n = 1;
for (LR::Basisfunction *b : lrspline->getElement(iel-1)->support() ) {
for (size_t i = 1; i <= 3; i++)
X(i,n) = b->cp(i-1);
n++;
}
int n = 1;
for (LR::Basisfunction* b : el->support())
X.fillColumn(n++,&(*b->cp()));
#if SP_DEBUG > 2
std::cout <<"\nCoordinates for element "<< iel << X << std::endl;
std::cout <<"\nCoordinates for element "<< iel << X << std::endl;
#endif
return true;
return true;
}
void ASMu3D::getNodalCoordinates (Matrix& X) const
{
const int n = lrspline->nBasisFunctions();
X.resize(3,n);
X.resize(3,lrspline->nBasisFunctions());
size_t inod = 1;
for(LR::Basisfunction *b : lrspline->getAllBasisfunctions() )
for (size_t i = 0; i < 3; i++)
X(i+1,inod++) = b->cp(i);
size_t inod = 1;
for (LR::Basisfunction* b : lrspline->getAllBasisfunctions())
X.fillColumn(inod++,&(*b->cp()));
}
Vec3 ASMu3D::getCoord (size_t inod) const
{
LR::Basisfunction* basis = lrspline->getBasisfunction(inod-1);
RealArray::const_iterator cit = basis->cp();
return Vec3(*cit,*(cit+1),*(cit+2));
}
bool ASMu3D::updateCoords (const Vector& displ)
{
std::cerr <<"ASMu3D::updateCoords not implemented properly yet" << std::endl;
exit(776654);
if (!lrspline) return true; // silently ignore empty patches
if (shareFE) return true;
if (displ.size() != 3*MLGN.size())
{
std::cerr <<" *** ASMu3D::updateCoords: Invalid dimension "
<< displ.size() <<" on displacement vector, should be "
<< 3*MLGN.size() << std::endl;
return false;
}
// lrspline->deform(displ,3);
return true;
}
bool ASMu3D::getOrder (int& p1, int& p2, int& p3) const
{
if (!lrspline) return false;
p1 = lrspline->order(0);
p2 = lrspline->order(1);
p3 = lrspline->order(2);
return true;
std::cerr <<" *** ASMu3D::updateCoords not implemented!"<< std::endl;
return false;
}
size_t ASMu3D::getNoBoundaryElms (char lIndex, char ldim) const
{
if (!lrspline) return 0;
/*TODO fix later (see ASMu2D::getNoBoundaryElms)
if (ldim < 1 && lIndex > 0)
return 1;
@@ -962,7 +782,7 @@ size_t ASMu3D::getNoBoundaryElms (char lIndex, char ldim) const
return n1*n2;
}
#endif
*/
return 0;
}
@@ -2099,144 +1919,7 @@ bool ASMu3D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
}
bool ASMu3D::evaluate (const ASMbase* input, const Vector& locVec, Vector& vec)
void ASMu3D::getBoundaryNodes (int lIndex, IntVec& nodes) const
{
std::cerr << "ASMu3D::evaluate(...) is not properly implemented yet :(" << std::endl;
exit(776654);
#if 0
ASMu3D* pch = (ASMu3D*)input;
// Compute parameter values of the result sampling points (Greville points)
RealArray gpar[3];
for (int dir = 0; dir < 3; dir++)
if (!this->getGrevilleParameters(gpar[dir],dir))
return false;
Matrix sValues;
pch->evalSolution(sValues, locVec, gpar, true);
// Project the results onto the spline basis to find control point
// values based on the result values evaluated at the Greville points.
// Note that we here implicitly assume that the number of Greville points
// equals the number of control points such that we don't have to resize
// the result array. Think that is always the case, but beware if trying
// other projection schemes later.
RealArray weights;
if (lrspline->rational())
lrspline->getWeights(weights);
const Vector& vec2 = sValues;
Go::SplineVolume* vol_new =
Go::VolumeInterpolator::regularInterpolation(lrspline->basis(0),
lrspline->basis(1),
lrspline->basis(2),
gpar[0], gpar[1], gpar[2],
const_cast<Vector&>(vec2),
sValues.rows(),
lrspline->rational(),
weights);
vec.resize(vol_new->coefs_end()-vol_new->coefs_begin());
std::copy(vol_new->coefs_begin(), vol_new->coefs_end(), vec.begin());
delete vol_new;
return true;
#endif
return false;
// TODO: Implement this before attempting FSI simulations with LR B-splines
}
bool ASMu3D::addXElms (short int dim, short int item, size_t nXn, IntVec& nodes)
{
std::cerr << "ASMu3D::addXElms not implemented properly yet\n";
exit(776654);
#if 0
if (dim != 2)
{
std::cerr <<" *** ASMu3D::addXElms: Invalid boundary dimension "<< dim
<<", only 2 (face) is allowed."<< std::endl;
return false;
}
else if (!svol || shareFE)
return false;
for (size_t i = 0; i < nXn; i++)
{
if (nodes.size() == i)
nodes.push_back(++gNod);
myMLGN.push_back(nodes[i]);
}
const int n1 = svol->numCoefs(0);
const int n2 = svol->numCoefs(1);
const int n3 = svol->numCoefs(2);
const int p1 = svol->order(0);
const int p2 = svol->order(1);
const int p3 = svol->order(2);
nXelm = (n1-p1+1)*(n2-p2+1)*(n3-p3+1);
myMNPC.resize(2*nXelm);
myMLGE.resize(2*nXelm,0);
int iel = 0;
bool skipMe = false;
for (int i3 = p3; i3 <= n3; i3++)
for (int i2 = p2; i2 <= n2; i2++)
for (int i1 = p1; i1 <= n1; i1++, iel++)
{
if (MLGE[iel] < 1) continue; // Skip zero-volume element
// Skip elements that are not on current boundary face
switch (item)
{
case 1: skipMe = i1 > p1; break;
case 2: skipMe = i1 < n1; break;
case 3: skipMe = i2 > p2; break;
case 4: skipMe = i2 < n2; break;
case 5: skipMe = i3 > p3; break;
case 6: skipMe = i3 < n3; break;
}
if (skipMe) continue;
IntVec& mnpc = myMNPC[nXelm+iel];
if (!mnpc.empty())
{
std::cerr <<" *** ASMu3D::addXElms: Only one X-face allowed."
<< std::endl;
return false;
}
mnpc = MNPC[iel]; // Copy the ordinary element nodes
// Negate node numbers that are not on the boundary face, to flag that
// they shall not receive any tangent and/or residual contributions
int lnod = 0;
for (int j3 = 0; j3 < p3; j3++)
for (int j2 = 0; j2 < p2; j2++)
for (int j1 = 0; j1 < p1; j1++, lnod++)
{
switch (item)
{
case 1: skipMe = j1 > 0; break;
case 2: skipMe = j1 < p1-1; break;
case 3: skipMe = j2 > 0; break;
case 4: skipMe = j2 < p2-1; break;
case 5: skipMe = j3 > 0; break;
case 6: skipMe = j3 < p3-1; break;
}
if (skipMe) // Hack for node 0: Using -maxint as flag instead
mnpc[lnod] = mnpc[lnod] == 0 ? -2147483648 : -mnpc[lnod];
}
// Add connectivity to the extra-ordinary nodes
for (size_t i = 0; i < nXn; i++)
mnpc.push_back(MLGN.size()-nXn+i);
myMLGE[nXelm+iel] = ++gEl;
}
#endif
return true;
}

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

@@ -14,447 +14,388 @@
#ifndef _ASM_U3D_H
#define _ASM_U3D_H
#include "ASM3D.h"
#include "ASMunstruct.h"
#include "ASM3D.h"
class FiniteElement;
namespace Go {
class SplineSurface;
class SplineVolume;
class SplineVolume;
}
namespace LR {
class LRSplineSurface;
class LRSplineVolume;
class LRSplineVolume;
}
/*!
\brief Driver for assembly of unstructured 3D spline FE models.
\details This class contains methods common for unstructured 3D spline patches.
\brief Driver for assembly of unstructured 3D spline FE models.
\details This class contains methods common for 3D LR-spline patches.
*/
class ASMu3D : public ASMunstruct, public ASM3D
{
public:
//! \brief Default constructor.
ASMu3D(unsigned char n_f = 3);
//! \brief Copy constructor.
ASMu3D(const ASMu3D& patch, unsigned char n_f = 0);
//! \brief Empty destructor.
virtual ~ASMu3D() {}
//! \brief Returns the spline volume representing the geometry of this patch.
LR::LRSplineVolume* getVolume() const { return lrspline; }
// Methods for model generation and refinement
// ===========================================
//! \brief Creates an instance by reading the given input stream.
virtual bool read(std::istream&);
//! \brief Writes the geometry of the SplineVolume object to given stream.
virtual bool write(std::ostream&, int = 0) const;
//! \brief Generates the finite element topology data for the patch.
//! \details The data generated are the element-to-node connectivity array,
//! and the arrays of global node and element numbers.
virtual bool generateFEMTopology();
//! \brief Clears the contents of the patch, making it empty.
//! \param[in] retainGeometry If \e true, the spline geometry is not cleared.
//! This is used to reinitialize the patch after it has been refined.
virtual void clear(bool retainGeometry = false);
//! \brief Returns a matrix with nodal coordinates for an element.
//! \param[in] iel Element index
//! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
//! in one element
virtual bool getElementCoordinates(Matrix& X, int iel) const;
//! \brief Returns a matrix with all nodal coordinates within the patch.
//! \param[out] X 3\f$\times\f$n-matrix, where \a n is the number of nodes
//! in the patch
virtual void getNodalCoordinates(Matrix& X) const;
//! \brief Returns the global coordinates for the given node.
//! \param[in] inod 1-based node index local to current patch
virtual Vec3 getCoord(size_t inod) const;
//! \brief Updates the nodal coordinates for this patch.
//! \param[in] displ Incremental displacements to update the coordinates with
virtual bool updateCoords(const Vector& displ);
//! \brief Finds the global numbers of the nodes on a patch boundary.
//! \param[in] lIndex Local index of the boundary edge
//! \param glbNodes Array of global boundary node numbers
virtual void getBoundaryNodes(int lIndex, IntVec& glbNodes) 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
virtual bool uniformRefine(int dir, int nInsert);
//! \brief Refines the parametrization by inserting extra tensor knots.
//! \param[in] dir Parameter direction to refine
//! \param[in] xi Relative positions of added knots in each existing knot span
virtual bool refine(int dir, const RealArray& xi);
//! \brief Raises the order of the tensor spline object for this patch.
//! \param[in] ru Number of times to raise the order in u-direction
//! \param[in] rv Number of times to raise the order in v-direction
//! \param[in] rw Number of times to raise the order in w-direction
virtual bool raiseOrder(int ru, int rv, int rw);
// Various methods for preprocessing of boundary conditions and patch topology
// ===========================================================================
//! \brief Constrains all DOFs on a given boundary face.
//! \param[in] dir Parameter direction defining the face to constrain
//! \param[in] open If \e true, exclude all points along the face boundary
//! \param[in] dof Which DOFs to constrain at each node on the face
//! \param[in] code Inhomogeneous dirichlet condition code
void constrainFace(int dir, bool open, int dof = 123, int code = 0);
//! \brief Constrains all DOFs in local directions on a given boundary face.
//! \param[in] dir Parameter direction defining the face to constrain
//! \param[in] open If \e true, exclude all points along the face boundary
//! \param[in] dof Which local DOFs to constrain at each node on the face
//! \param[in] code Inhomogeneous dirichlet condition code
//! \param[in] project If \e true, the local axis directions are projected
//! \param[in] T1 Desired global direction of first local tangent direction
//! \return Number of additional nodes added due to local axis constraints
size_t constrainFaceLocal(int dir, bool open, int dof = 3, int code = 0,
bool project = false, char T1 = '\0');
//! \brief Constrains all DOFs on a given boundary edge.
//! \param[in] lEdge Local index [1,12] of the edge to constrain
//! \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
void constrainEdge(int lEdge, bool open, int dof = 123, int code = 0);
//! \brief Constrains all DOFs along a line on a given boundary face.
//! \param[in] fdir Parameter direction defining the face to constrain
//! \param[in] ldir Parameter direction defining the line to constrain
//! \param[in] xi Parameter value defining the line to constrain
//! \param[in] dof Which DOFs to constrain at each node along the line
//! \param[in] code Inhomogeneous dirichlet condition code
//!
//! \details The parameter \a xi has to be in the domain [0.0,1.0], where
//! 0.0 means the beginning of the domain and 1.0 means the end. The line to
//! constrain goes along the parameter direction \a ldir in the face with
//! with normal in parameter direction \a fdir, and positioned along the third
//! parameter direction as indicated by \a xi. The actual value of \a xi
//! is converted to the integer value closest to \a xi*n, where \a n is the
//! number of nodes (control points) in that parameter direction.
void constrainLine(int fdir, int ldir, double xi,
int dof = 123, int code = 0);
//! \brief Constrains a corner node identified by the three parameter indices.
//! \param[in] I Parameter index in u-direction
//! \param[in] J Parameter index in v-direction
//! \param[in] K Parameter index in w-direction
//! \param[in] dof Which DOFs to constrain at the node
//! \param[in] code Inhomogeneous dirichlet condition code
//!
//! \details The sign of the three indices is used to define whether we want
//! the node at the beginning or the end of that parameter direction.
//! The magnitude of the indices are not used.
void constrainCorner(int I, int J, int K,
int dof = 123, int code = 0);
//! \brief Constrains a node identified by three relative parameter values.
//! \param[in] xi Parameter in u-direction
//! \param[in] eta Parameter in v-direction
//! \param[in] zeta Parameter in w-direction
//! \param[in] dof Which DOFs to constrain at the node
//! \param[in] code Inhomogeneous dirichlet condition code
//!
//! \details The parameter values have to be in the domain [0.0,1.0], where
//! 0.0 means the beginning of the domain and 1.0 means the end. For values
//! in between, the actual index is taken as the integer value closest to
//! \a r*n, where \a r denotes the given relative parameter value,
//! and \a n is the number of nodes along that parameter direction.
void constrainNode(double xi, double eta, double zeta,
int dof = 123, int code = 0);
/* More multipatch stuff, maybe later...
//! \brief Connects all matching nodes on two adjacent boundary faces.
//! \param[in] face Local face index of this patch, in range [1,6]
//! \param neighbor The neighbor patch
//! \param[in] nface Local face index of neighbor patch, in range [1,6]
//! \param[in] norient Relative face orientation flag (see below)
//!
//! \details The face orientation flag \a norient must be in range [0,7].
//! When interpreted as a binary number, its 3 digits are decoded as follows:
//! - left digit = 1: The u and v parameters of the neighbor face are swapped
//! - middle digit = 1: Parameter \a u in neighbor patch face is reversed
//! - right digit = 1: Parameter \a v in neighbor patch face is reversed
virtual bool connectPatch(int face, ASMu3D& neighbor, int nface,
int norient = 0);
//! \brief Makes two opposite boundary faces periodic.
//! \param[in] dir Parameter direction defining the periodic faces
//! \param[in] basis Which basis to connect (mixed methods), 0 means both
//! \param[in] master 1-based index of the first master node in this basis
virtual void closeFaces(int dir, int basis = 0, int master = 1);
*/
// Methods for integration of finite element quantities.
// These are the main computational methods of the ASM class hierarchy.
// ====================================================================
//! \brief Evaluates an integral over the interior patch domain.
//! \param integrand Object with problem-specific data and methods
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
virtual bool integrate(Integrand& integrand,
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch face.
//! \param integrand Object with problem-specific data and methods
//! \param[in] lIndex Local index [1,6] of the boundary face
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
virtual bool integrate(Integrand& integrand, int lIndex,
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch edge.
//! \param integrand Object with problem-specific data and methods
//! \param[in] lEdge Local index [1,12] of the patch edge
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
virtual bool integrateEdge(Integrand& integrand, int lEdge,
GlobalIntegral& glbInt, const TimeDomain& time);
// Post-processing methods
// =======================
//! \brief Evaluates the geometry at a specified point.
//! \param[in] xi Dimensionless parameters in range [0.0,1.0] of the point
//! \param[out] param The (u,v,w) parameters of the point in knot-span domain
//! \param[out] X The Cartesian coordinates of the point
//! \return Local node number within the patch that matches the point, if any
//! \return 0 if no node (control point) matches this point
virtual int evalPoint(const double* xi, double* param, Vec3& X) const;
//! \brief Calculates parameter values for visualization nodal points.
//! \param[out] prm Parameter values in given direction for all points
//! \param[in] dir Parameter direction (0,1,2)
//! \param[in] nSegSpan Number of visualization segments over each knot-span
virtual bool getGridParameters(RealArray& prm, int dir, int nSegSpan) const;
//! \brief Creates a hexahedron element model of this patch for visualization.
//! \param[out] grid The generated hexahedron grid
//! \param[in] npe Number of visualization nodes over each knot span
//! \note The number of element nodes must be set in \a grid on input.
virtual bool tesselate(ElementBlock& grid, const int* npe) const;
//! \brief Evaluates the primary solution field at all visualization points.
//! \param[out] sField Solution field
//! \param[in] locSol Solution vector in DOF-order
//! \param[in] npe Number of visualization nodes over each knot span
virtual bool evalSolution(Matrix& sField, const Vector& locSol,
const int* npe) const;
//! \brief Evaluates the primary solution field at the given points.
//! \param[out] sField Solution field
//! \param[in] locSol Solution vector local to current patch
//! \param[in] gpar Parameter values of the result sampling points
//! \param[in] regular Flag indicating how the sampling points are defined
//!
//! \details When \a regular is \e true, it is assumed that the parameter
//! value array \a gpar forms a regular tensor-product point grid of dimension
//! \a gpar[0].size() \a X \a gpar[1].size() \a X \a gpar[2].size().
//! Otherwise, we assume that it contains the \a u, \a v and \a w parameters
//! directly for each sampling point.
virtual bool evalSolution(Matrix& sField, const Vector& locSol,
const RealArray* gpar, bool regular = true) const;
//! \brief Evaluates the secondary solution field at all visualization points.
//! \param[out] sField Solution field
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] npe Number of visualization nodes over each knot span
//! \param[in] project Flag indicating result recovery method
//! (0=none, 'D'=direct evaluation, 'S'=superconvergent recovery)
//!
//! \details The secondary solution is derived from the primary solution,
//! which is assumed to be stored within the \a integrand for current patch.
//! If \a npe is NULL, the solution is recovered or evaluated at the Greville
//! points and then projected onto the spline basis to obtain the control
//! point values, which then are returned through \a sField.
//! If \a npe is not NULL and \a project is defined, the solution is also
//! projected onto the spline basis, and then evaluated at the \a npe points.
virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand,
const int* npe = 0, char project = '\0') const;
public:
//! \brief Default constructor.
ASMu3D(unsigned char n_f = 3);
//! \brief Copy constructor.
ASMu3D(const ASMu3D& patch, unsigned char n_f = 0);
//! \brief Empty destructor.
virtual ~ASMu3D() {}
//! \brief Projects the secondary solution field onto the primary basis.
//! \param[in] integrand Object with problem-specific data and methods
virtual LR::LRSpline* evalSolution(const IntegrandBase& integrand) const { return NULL; }
//! \brief Returns the spline volume representing the geometry of this patch.
LR::LRSplineVolume* getVolume() const { return lrspline; }
//! \brief Returns the spline surface representing a boundary of this patch.
//! \param[in] dir Parameter direction defining which boundary to return
virtual LR::LRSplineSurface* getBoundary(int dir);
//! \brief Returns the spline volume representing the basis of this patch.
virtual LR::LRSplineVolume* getBasis(int = 1) const { return lrspline; }
//! \brief Evaluates the secondary solution field at the given points.
//! \param[out] sField Solution field
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] gpar Parameter values of the result sampling points
//! \param[in] regular Flag indicating how the sampling points are defined
//!
//! \details The secondary solution is derived from the primary solution,
//! which is assumed to be stored within the \a integrand for current patch.
//! When \a regular is \e true, it is assumed that the parameter value array
//! \a gpar forms a regular tensor-product point grid of dimension
//! \a gpar[0].size() \a X \a gpar[1].size() \a X \a gpar[2].size().
//! Otherwise, we assume that it contains the \a u, \a v and \a w parameters
//! directly for each sampling point.
virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand,
const RealArray* gpar, bool regular = true) const;
// Methods for model generation
// ============================
//! \brief Creates an instance by reading the given input stream.
virtual bool read(std::istream&);
//! \brief Writes the geometry of the SplineVolume object to given stream.
virtual bool write(std::ostream&, int = 0) const;
//! \brief Generates the finite element topology data for the patch.
//! \details The data generated are the element-to-node connectivity array,
//! the node-to-IJK-index array, as well as global node and element numbers.
virtual bool generateFEMTopology();
//! \brief Clears the contents of the patch, making it empty.
//! \param[in] retainGeometry If \e true, the spline geometry is not cleared.
//! This is used to reinitialize the patch after it has been refined.
virtual void clear(bool retainGeometry = false);
//! \brief Adds extraordinary elements associated with a patch boundary.
//! \param[in] dim Dimension of the boundary (should be 2)
//! \param[in] item Local index of the boundary face
//! \param[in] nXn Number of extraordinary nodes
//! \param[out] nodes Global numbers assigned to the extraordinary nodes
virtual bool addXElms(short int dim, short int item,
size_t nXn, IntVec& nodes);
//! \brief Returns local 1-based index of the node with given global number.
//! \details If the given node number is not present, 0 is returned.
//! \param[in] globalNum Global node number
//! \param[in] noAddedNodes If \e true, use \a xnMap to find the real node
virtual size_t getNodeIndex(int globalNum, bool noAddedNodes = false) const;
//! \brief Returns the total number of nodes in this patch.
virtual size_t getNoNodes(int basis = 0) const;
//! \brief Returns a matrix with nodal coordinates for an element.
//! \param[in] iel Element index
//! \param[out] X 3\f$\times\f$n-matrix,
//! where \a n is the number of nodes in one element
virtual bool getElementCoordinates(Matrix& X, int iel) const;
//! \brief Returns a matrix with all nodal coordinates within the patch.
//! \param[out] X 3\f$\times\f$n-matrix,
//! where \a n is the number of nodes in the patch
virtual void getNodalCoordinates(Matrix& X) const;
//! \brief Returns the global coordinates for the given node.
//! \param[in] inod 1-based node index local to current patch
virtual Vec3 getCoord(size_t inod) const;
//! \brief Updates the nodal coordinates for this patch.
//! \param[in] displ Incremental displacements to update the coordinates with
virtual bool updateCoords(const Vector& displ);
//! \brief Checks that the patch is modelled in a right-hand-side system.
virtual bool checkRightHandSystem();
//! \brief Refines the parametrization by inserting extra knots.
//! \param[in] dir Parameter direction to refine
//! \param[in] xi Relative positions of added knots in each existing knot span
bool refine(int dir, const RealArray& xi);
//! \brief Refines the parametrization by inserting extra knots uniformly.
//! \param[in] dir Parameter direction to refine
//! \param[in] nInsert Number of extra knots to insert in each knot-span
bool uniformRefine(int dir, int nInsert);
//! \brief Raises the order of the SplineVolume object for this patch.
//! \param[in] ru Number of times to raise the order in u-direction
//! \param[in] rv Number of times to raise the order in v-direction
//! \param[in] rw Number of times to raise the order in w-direction
bool raiseOrder(int ru, int rv, int rw);
bool refine(const std::vector<int>& elements,
const std::vector<int>& options, const char* fName = 0) {
std::cout << "ASMu3D::refine()" << std::endl;
return ASMunstruct::refine(elements, options, fName);
}
// Various methods for preprocessing of boundary conditions and patch topology
// ===========================================================================
//! \brief Constrains all DOFs on a given boundary face.
//! \param[in] dir Parameter direction defining the face to constrain
//! \param[in] open If \e true, exclude all points along the face boundary
//! \param[in] dof Which DOFs to constrain at each node on the face
//! \param[in] code Inhomogeneous dirichlet condition code
void constrainFace(int dir, bool open, int dof = 123, int code = 0);
//! \brief Constrains all DOFs in local directions on a given boundary face.
//! \param[in] dir Parameter direction defining the face to constrain
//! \param[in] open If \e true, exclude all points along the face boundary
//! \param[in] dof Which local DOFs to constrain at each node on the face
//! \param[in] code Inhomogeneous dirichlet condition code
//! \param[in] project If \e true, the local axis directions are projected
//! \param[in] T1 Desired global direction of first local tangent direction
//! \return Number of additional nodes added due to local axis constraints
size_t constrainFaceLocal(int dir, bool open, int dof = 3, int code = 0,
bool project = false, char T1 = '\0');
//! \brief Constrains all DOFs on a given boundary edge.
//! \param[in] lEdge Local index [1,12] of the edge to constrain
//! \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
void constrainEdge(int lEdge, bool open, int dof = 123, int code = 0);
//! \brief Constrains all DOFs along a line on a given boundary face.
//! \param[in] fdir Parameter direction defining the face to constrain
//! \param[in] ldir Parameter direction defining the line to constrain
//! \param[in] xi Parameter value defining the line to constrain
//! \param[in] dof Which DOFs to constrain at each node along the line
//! \param[in] code Inhomogeneous dirichlet condition code
//!
//! \details The parameter \a xi has to be in the domain [0.0,1.0], where
//! 0.0 means the beginning of the domain and 1.0 means the end. The line to
//! constrain goes along the parameter direction \a ldir in the face with
//! with normal in parameter direction \a fdir, and positioned along the third
//! parameter direction as indicated by \a xi. The actual value of \a xi
//! is converted to the integer value closest to \a xi*n, where \a n is the
//! number of nodes (control points) in that parameter direction.
void constrainLine(int fdir, int ldir, double xi,
int dof = 123, int code = 0);
//! \brief Constrains a corner node identified by the three parameter indices.
//! \param[in] I Parameter index in u-direction
//! \param[in] J Parameter index in v-direction
//! \param[in] K Parameter index in w-direction
//! \param[in] dof Which DOFs to constrain at the node
//! \param[in] code Inhomogeneous dirichlet condition code
//!
//! \details The sign of the three indices is used to define whether we want
//! the node at the beginning or the end of that parameter direction.
//! The magnitude of the indices are not used.
void constrainCorner(int I, int J, int K,
int dof = 123, int code = 0);
//! \brief Constrains a node identified by three relative parameter values.
//! \param[in] xi Parameter in u-direction
//! \param[in] eta Parameter in v-direction
//! \param[in] zeta Parameter in w-direction
//! \param[in] dof Which DOFs to constrain at the node
//! \param[in] code Inhomogeneous dirichlet condition code
//!
//! \details The parameter values have to be in the domain [0.0,1.0], where
//! 0.0 means the beginning of the domain and 1.0 means the end. For values
//! in between, the actual index is taken as the integer value closest to
//! \a r*n, where \a r denotes the given relative parameter value,
//! and \a n is the number of nodes along that parameter direction.
void constrainNode(double xi, double eta, double zeta,
int dof = 123, int code = 0);
//! \brief Connects all matching nodes on two adjacent boundary faces.
//! \param[in] face Local face index of this patch, in range [1,6]
//! \param neighbor The neighbor patch
//! \param[in] nface Local face index of neighbor patch, in range [1,6]
//! \param[in] norient Relative face orientation flag (see below)
//!
//! \details The face orientation flag \a norient must be in range [0,7].
//! When interpreted as a binary number, its 3 digits are decoded as follows:
//! - left digit = 1: The u and v parameters of the neighbor face are swapped
//! - middle digit = 1: Parameter \a u in neighbor patch face is reversed
//! - right digit = 1: Parameter \a v in neighbor patch face is reversed
virtual bool connectPatch(int face, ASMu3D& neighbor, int nface,
int norient = 0);
//! \brief Makes two opposite boundary faces periodic.
//! \param[in] dir Parameter direction defining the periodic faces
//! \param[in] basis Which basis to connect (mixed methods), 0 means both
//! \param[in] master 1-based index of the first master node in this basis
virtual void closeFaces(int dir, int basis = 0, int master = 1);
//! \brief Updates the time-dependent in-homogeneous Dirichlet coefficients.
//! \param[in] func Scalar property fields
//! \param[in] vfunc Vector property fields
//! \param[in] time Current time
virtual bool updateDirichlet(const std::map<int,RealFunc*>& func,
const std::map<int,VecFunc*>& vfunc,
double time = 0.0);
// Methods for integration of finite element quantities.
// These are the main computational methods of the ASM class hierarchy.
// ====================================================================
//! \brief Evaluates an integral over the interior patch domain.
//! \param integrand Object with problem-specific data and methods
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
virtual bool integrate(Integrand& integrand,
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch face.
//! \param integrand Object with problem-specific data and methods
//! \param[in] lIndex Local index [1,6] of the boundary face
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
virtual bool integrate(Integrand& integrand, int lIndex,
GlobalIntegral& glbInt, const TimeDomain& time);
//! \brief Evaluates a boundary integral over a patch edge.
//! \param integrand Object with problem-specific data and methods
//! \param[in] lEdge Local index [1,12] of the patch edge
//! \param glbInt The integrated quantity
//! \param[in] time Parameters for nonlinear/time-dependent simulations
virtual bool integrateEdge(Integrand& integrand, int lEdge,
GlobalIntegral& glbInt, const TimeDomain& time);
// Post-processing methods
// =======================
//! \brief Evaluates the geometry at a specified point.
//! \param[in] xi Dimensionless parameters in range [0.0,1.0] of the point
//! \param[out] param The (u,v,w) parameters of the point in knot-span domain
//! \param[out] X The Cartesian coordinates of the point
//! \return Local node number within the patch that matches the point, if any
//! \return 0 if no node (control point) matches this point
virtual int evalPoint(const double* xi, double* param, Vec3& X) const;
//! \brief Calculates parameter values for visualization nodal points.
//! \param[out] prm Parameter values in given direction for all points
//! \param[in] dir Parameter direction (0,1,2)
//! \param[in] nSegSpan Number of visualization segments over each knot-span
virtual bool getGridParameters(RealArray& prm, int dir, int nSegSpan) const;
//! \brief Creates a hexahedron element model of this patch for visualization.
//! \param[out] grid The generated hexahedron grid
//! \param[in] npe Number of visualization nodes over each knot span
//! \note The number of element nodes must be set in \a grid on input.
virtual bool tesselate(ElementBlock& grid, const int* npe) const;
//! \brief Evaluates the primary solution field at all visualization points.
//! \param[out] sField Solution field
//! \param[in] locSol Solution vector in DOF-order
//! \param[in] npe Number of visualization nodes over each knot span
virtual bool evalSolution(Matrix& sField, const Vector& locSol,
const int* npe) const;
//! \brief Evaluates the primary solution field at the given points.
//! \param[out] sField Solution field
//! \param[in] locSol Solution vector local to current patch
//! \param[in] gpar Parameter values of the result sampling points
//! \param[in] regular Flag indicating how the sampling points are defined
//!
//! \details When \a regular is \e true, it is assumed that the parameter
//! value array \a gpar forms a regular tensor-product point grid of dimension
//! \a gpar[0].size() \a X \a gpar[1].size() \a X \a gpar[2].size().
//! Otherwise, we assume that it contains the \a u, \a v and \a w parameters
//! directly for each sampling point.
virtual bool evalSolution(Matrix& sField, const Vector& locSol,
const RealArray* gpar, bool regular = true) const;
//! \brief Evaluates the secondary solution field at all visualization points.
//! \param[out] sField Solution field
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] npe Number of visualization nodes over each knot span
//! \param[in] project Flag indicating result recovery method
//! (0=none, 'D'=direct evaluation, 'S'=superconvergent recovery)
//!
//! \details The secondary solution is derived from the primary solution,
//! which is assumed to be stored within the \a integrand for current patch.
//! If \a npe is NULL, the solution is recovered or evaluated at the Greville
//! points and then projected onto the spline basis to obtain the control
//! point values, which then are returned through \a sField.
//! If \a npe is not NULL and \a project is defined, the solution is also
//! projected onto the spline basis, and then evaluated at the \a npe points.
virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand,
const int* npe = 0, char project = '\0') const;
public:
//! \brief Projects the secondary solution field onto the primary basis.
//! \param[in] integrand Object with problem-specific data and methods
virtual LR::LRSpline* evalSolution(const IntegrandBase& integrand) const { return NULL; }
//! \brief Evaluates the secondary solution field at the given points.
//! \param[out] sField Solution field
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] gpar Parameter values of the result sampling points
//! \param[in] regular Flag indicating how the sampling points are defined
//!
//! \details The secondary solution is derived from the primary solution,
//! which is assumed to be stored within the \a integrand for current patch.
//! When \a regular is \e true, it is assumed that the parameter value array
//! \a gpar forms a regular tensor-product point grid of dimension
//! \a gpar[0].size() \a X \a gpar[1].size() \a X \a gpar[2].size().
//! Otherwise, we assume that it contains the \a u, \a v and \a w parameters
//! directly for each sampling point.
virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand,
const RealArray* gpar, bool regular = true) const;
//! \brief Projects the secondary solution using a discrete global L2-norm.
//! \param[out] sField Secondary solution field control point values
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] continuous If \e true, a continuous L2-projection is used
virtual bool globalL2projection(Matrix& sField,
const IntegrandBase& integrand,
bool continuous = false) const { return false; }
//! \brief Evaluate and interpolate a field over a given geometry
//! \param[in] input The basis of the field to evaluate
//! \param[in] locVec The coefficients of the field
//! \param[out] vec The obtained coefficients after interpolation
virtual bool evaluate(const ASMbase* input,
const Vector& locVec, Vector& vec);
//! \brief Projects the secondary solution using a discrete global L2-norm.
//! \param[out] sField Secondary solution field control point values
//! \param[in] integrand Object with problem-specific data and methods
//! \param[in] continuous If \e true, a continuous L2-projection is used
virtual bool globalL2projection(Matrix& sField,
const IntegrandBase& integrand,
bool continuous = false) const { return false; }
protected:
// Internal utility methods
// ========================
// Internal utility methods
// ========================
//! \brief Connects all matching nodes on two adjacent boundary faces.
//! \param[in] face Local face index of this patch, in range [1,6]
//! \param neighbor The neighbor patch
//! \param[in] nface Local face index of neighbor patch, in range [1,6]
//! \param[in] norient Relative face orientation flag (see \a connectPatch)
//! \param[in] basis Which basis to connect the nodes for (mixed methods)
//! \param[in] slave 0-based index of the first slave node in this basis
//! \param[in] master 0-based index of the first master node in this basis
bool connectBasis(int face, ASMu3D& neighbor, int nface, int norient,
int basis = 1, int slave = 0, int master = 0);
//! \brief Connects all matching nodes on two adjacent boundary faces.
//! \param[in] face Local face index of this patch, in range [1,6]
//! \param neighbor The neighbor patch
//! \param[in] nface Local face index of neighbor patch, in range [1,6]
//! \param[in] norient Relative face orientation flag (see \a connectPatch)
//! \param[in] basis Which basis to connect the nodes for (mixed methods)
//! \param[in] slave 0-based index of the first slave node in this basis
//! \param[in] master 0-based index of the first master node in this basis
bool connectBasis(int face, ASMu3D& neighbor, int nface, int norient,
int basis = 1, int slave = 0, int master = 0);
//! \brief Extracts parameter values of the Gauss points in one direction.
//! \param[out] uGP Parameter values in given direction for all points
//! \param[in] dir Parameter direction (0,1,2)
//! \param[in] nGauss Number of Gauss points along a knot-span
//! \param[in] iel Element index
//! \param[in] xi Dimensionless Gauss point coordinates [-1,1]
//! \return The parameter value matrix casted into a one-dimensional vector
void getGaussPointParameters(RealArray& uGP, int dir, int nGauss,
int iel, const double* xi) const;
//! \brief Extracts parameter values of the Gauss points in one direction.
//! \param[out] uGP Parameter values in given direction for all points
//! \param[in] dir Parameter direction (0,1,2)
//! \param[in] nGauss Number of Gauss points along a knot-span
//! \param[in] iel Element index
//! \param[in] xi Dimensionless Gauss point coordinates [-1,1]
//! \return The parameter value matrix casted into a one-dimensional vector
void getGaussPointParameters(RealArray& uGP, int dir, int nGauss,
int iel, const double* xi) const;
//! \brief Calculates parameter values for the Greville points.
//! \param[out] prm Parameter values in given direction for all points
//! \param[in] dir Parameter direction (0,1,2)
bool getGrevilleParameters(RealArray& prm, int dir) const;
//! \brief Calculates parameter values for the Greville points.
//! \param[out] prm Parameter values in given direction for all points
//! \param[in] dir Parameter direction (0,1,2)
bool getGrevilleParameters(RealArray& prm, int dir) const;
//! \brief Calculates parameter values for the Quasi-Interpolation points.
//! \param[out] prm Parameter values in given direction for all points
//! \param[in] dir Parameter direction (0,1,2)
bool getQuasiInterplParameters(RealArray& prm, int dir) const;
//! \brief Calculates parameter values for the Quasi-Interpolation points.
//! \param[out] prm Parameter values in given direction for all points
//! \param[in] dir Parameter direction (0,1,2)
bool getQuasiInterplParameters(RealArray& prm, int dir) const;
//! \brief Returns the volume in the parameter space for an element.
//! \param[in] iel Element index
double getParametricVolume(int iel) const;
//! \brief Returns the volume in the parameter space for an element.
//! \param[in] iel Element index
double getParametricVolume(int iel) const;
//! \brief Returns boundary face area in the parameter space for an element.
//! \param[in] iel Element index
//! \param[in] dir Local face index of the boundary face
double getParametricArea(int iel, int dir) const;
//! \brief Returns boundary face area in the parameter space for an element.
//! \param[in] iel Element index
//! \param[in] dir Local face index of the boundary face
double getParametricArea(int iel, int dir) const;
//! \brief Computes the element corner coordinates.
//! \param[in] iel index (0 based) of the element on which corners are requested
//! \param[out] XC Coordinates of the element corners
void getElementCorners(int iel, std::vector<Vec3>& XC) const;
//! \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;
//! \brief Evaluate all basis functions and \a derivs number of derivatives on one element
virtual void evaluateBasis(FiniteElement &el, int derivs) const;
//! \brief Evaluate all basis functions and \a derivs number of derivatives on one element
virtual void evaluateBasis(FiniteElement &el, int derivs) const;
//! \brief Evaluate all basis functions and first derivatives on one element
virtual void evaluateBasis(FiniteElement &el, Matrix &dNdu, Matrix &C, Matrix &B) const ;
//! \brief Evaluate all basis functions and first derivatives on one element
virtual void evaluateBasis(FiniteElement &el, Matrix &dNdu, Matrix &C, Matrix &B) const ;
//! \brief Evaluate all basis functions and first derivatives on one element
virtual void evaluateBasis(FiniteElement &el, Matrix &dNdu) const;
//! \brief Evaluate all basis functions and first derivatives on one element
virtual void evaluateBasis(FiniteElement &el, Matrix &dNdu) const;
//! \brief Evaluate all basis functions and second order derivatives on one element
virtual void evaluateBasis(FiniteElement &el, Matrix &dNdu, Matrix3D& d2Ndu2) const;
//! \brief Evaluate all basis functions and second order derivatives on one element
virtual void evaluateBasis(FiniteElement &el, Matrix &dNdu, Matrix3D& d2Ndu2) const;
public:
//! \brief Returns the polynomial order in each parameter direction.
//! \param[out] p1 Order in first (u) direction
//! \param[out] p2 Order in second (v) direction
//! \param[out] p3 Order in third (w) direction
virtual bool getOrder(int& p1, int& p2, int& p3) const;
//! \brief Returns the number of elements on a boundary.
virtual size_t getNoBoundaryElms(char lIndex, char ldim) const;
// int getNodeID(size_t inod, bool = false) const { return inod; };
// unsigned char getNodalDOFs(size_t inod) const { return 0; };
// void getNoIntPoints(size_t& nPt) { };
// void getNoBouPoints(size_t& nPt, char ldim, char lindx) { };
// bool getSolution(Matrix& sField, const Vector& locSol,
// const IntVec& nodes) const { return false; };
// void extractNodeVec(const Vector& globVec, Vector& nodeVec,
// unsigned char nndof = 0, int basis = 1) const { };
// bool injectNodeVec(const Vector& nodeVec, Vector& globVec,
// unsigned char nndof = 0) const { return false;};
//! \brief Returns the number of elements on a boundary.
virtual size_t getNoBoundaryElms(char lIndex, char ldim) const;
protected:
LR::LRSplineVolume* lrspline; //!< Pointer to the actual spline volume object
LR::LRSplineVolume* lrspline; //!< Pointer to the LR-spline volume object
std::map<size_t,size_t> xnMap; //!< Node index map used by \a getCoord
std::map<size_t,size_t> nxMap; //!< Node index map used by \a getNodeID
std::vector<Matrix> bezierExtract; //!< The Bezier exctraction matrices for all elements
std::vector<Matrix> bezierExtract; //!< Bezier extraction matrices
private:
Go::SplineVolume* tensorspline; //!< Pointer to original tensor spline object
// The tensor spline object is kept for backward compatability with the REFINE
// and RAISEORDER key-words, although we take note that there is a possibility
// of optimization since all mapping values and Jacobians may be performed on
// this object for increased efficiency.
Go::SplineVolume* tensorspline; //!< Pointer to original tensor spline object
// The tensor spline object is kept for backward compatability with the REFINE
// and RAISEORDER key-words, although we take note that there is a possibility
// of optimization since all mapping values and Jacobians may be performed on
// this object for increased efficiency.
};
#endif

83
src/SIM/ForceIntegrator.C
View File

@@ -7,7 +7,7 @@
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Driver for integration of boundary forces.
//! \brief Driver for integration of boundary and nodal forces.
//!
//==============================================================================
@@ -15,7 +15,7 @@
#include "SIMbase.h"
#include "ASMbase.h"
#include "IntegrandBase.h"
#include "GlobalIntegral.h"
#include "GlbForceVec.h"
#ifdef PARALLEL_PETSC
#include "petscversion.h"
#if PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 2
@@ -28,7 +28,19 @@
Vector SIM::getBoundaryForce (const Vectors& solution, SIMbase* model, int code,
const TimeDomain& time, const Vec3* X0)
const TimeDomain& time, const Vector* X0)
{
if (X0)
{
Vec3 X(X0->ptr(),X0->size());
return SIM::getBoundaryForce(solution,model,code,time,&X);
}
return SIM::getBoundaryForce(solution,model,code,time);
}
Vector SIM::getBoundaryForce (const Vectors& solution, SIMbase* model, int code,
const TimeDomain& time, const Vec3* X0)
{
ForceBase* forceInt = model->getBoundaryForceIntegrand(X0);
if (!forceInt)
@@ -39,7 +51,7 @@ Vector SIM::getBoundaryForce (const Vectors& solution, SIMbase* model, int code,
forceInt->initBuffer(model->getNoElms());
const std::vector<ASMbase*>& feModel = model->getFEModel();
const ASMVec& feModel = model->getFEModel();
// Integrate forces for given boundary segment
bool ok = true;
@@ -77,6 +89,67 @@ Vector SIM::getBoundaryForce (const Vectors& solution, SIMbase* model, int code,
if (ok) return force;
std::cerr <<" *** SIM::getBoundaryForce: Failed to evaluate boundary force"
<< std::endl;
<< std::endl;
return Vector();
}
bool SIM::getNodalForces (const Vectors& solution, SIMbase* model, int code,
const TimeDomain& time, GlbForceVec& force)
{
force.initialize();
ForceBase* forceInt = model->getNodalForceIntegrand();
if (!forceInt)
{
std::cerr <<" *** SIM::getNodalForces: No force integrand."<< std::endl;
return false;
}
const ASMVec& feModel = model->getFEModel();
// Integrate nodal forces for given boundary segment
bool ok = true;
size_t prevPatch = 0;
PropertyVec::const_iterator p;
for (p = model->begin_prop(); p != model->end_prop() && ok; p++)
if (abs(p->pindx) == code)
{
size_t j = p->patch;
if (j < 1 || j > feModel.size())
ok = false;
else if (abs(p->ldim)+1 == feModel[j-1]->getNoSpaceDim())
{
if (j != prevPatch)
ok = model->extractPatchSolution(solution,j-1);
ok &= feModel[j-1]->integrate(*forceInt,abs(p->lindx),force,time);
prevPatch = j;
}
}
delete forceInt;
if (ok & force.finalize())
return true;
std::cerr <<" *** SIM::getNodalForces: Failed to evaluate nodal forces"
<< std::endl;
return false;
}
bool SIM::initBoundaryNodeMap (SIMbase* model, int code, GlbForceVec& force)
{
const ASMVec& feModel = model->getFEModel();
IntVec glbNodes;
PropertyVec::const_iterator p;
for (p = model->begin_prop(); p != model->end_prop(); p++)
if (abs(p->pindx) == code && p->patch > 0 && p->patch < feModel.size())
{
ASMbase* patch = feModel[p->patch-1];
if (abs(p->ldim)+1 == patch->getNoSpaceDim())
patch->getBoundaryNodes(abs(p->lindx),glbNodes);
}
return force.initNodeMap(glbNodes,model->getNoSpaceDim());
}

32
src/SIM/ForceIntegrator.h
View File

@@ -7,7 +7,7 @@
//!
//! \author Knut Morten Okstad / SINTEF
//!
//! \brief Driver for integration of boundary forces.
//! \brief Driver for integration of boundary and nodal forces.
//!
//==============================================================================
@@ -19,18 +19,44 @@
class SIMbase;
class TimeDomain;
class Vec3;
class GlbForceVec;
namespace SIM
{
//! \brief Integrate the force resultant on a specified boundary.
//! \brief Integrates the force resultant on a specified boundary.
//! \param[in] solution Primary solution vectors in DOF order
//! \param[in] model The isogeometric finite element model
//! \param[in] code Code indentifying the boundary subjected to integration
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] X0 Pivot point for torque calculation
//! \return The force (and torque) resultant
Vector getBoundaryForce(const Vectors& solution, SIMbase* model, int code,
const TimeDomain& time, const Vec3* X0 = NULL);
const TimeDomain& time, const Vector* X0);
//! \brief Integrates the force resultant on a specified boundary.
//! \param[in] solution Primary solution vectors in DOF order
//! \param[in] model The isogeometric finite element model
//! \param[in] code Code indentifying the boundary subjected to integration
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param[in] X0 Pivot point for torque calculation
//! \return The force (and torque) resultant
Vector getBoundaryForce(const Vectors& solution, SIMbase* model, int code,
const TimeDomain& time, const Vec3* X0 = NULL);
//! \brief Integrates nodal forces on a specified boundary.
//! \param[in] solution Primary solution vectors in DOF order
//! \param[in] model The isogeometric finite element model
//! \param[in] code Code indentifying the boundary subjected to integration
//! \param[in] time Parameters for nonlinear and time-dependent simulations
//! \param force Global nodal force container (compressed storage)
bool getNodalForces(const Vectors& solution, SIMbase* model, int code,
const TimeDomain& time, GlbForceVec& force);
//! \brief Detects the global nodes that reside on a specified boundary.
//! \param[in] model The isogeometric finite element model
//! \param[in] code Property code associated with the boundary
//! \param[out] force Global nodal force container (compressed storage)
bool initBoundaryNodeMap(SIMbase* model, int code, GlbForceVec& force);
}
#endif