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

Changed: Reimplemented SAMpatchPara::initConstraintEqs (which I suspect did

Browse Source 
not work anyway, not tested on examples with multi-point constraints) using
the parent class method and then expanding the mpmceq array to dof-size.

Added: virtual method ASMstruct::getSize(int&,int&,int&,int) and using that
in all the getSubdomains methods in SAMpatchPara (further casting not needed).

git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@2355 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo
2013-05-02 13:16:44 +00:00
committed by Knut Morten Okstad
parent 8c0579dcdc
commit 8c0f8fbe6f
9 changed files with 169 additions and 405 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
src/ASM/ASMs1D.C
View File

@@ -462,6 +462,14 @@ bool ASMs1D::getOrder (int& p1, int& p2, int& p3) const
}
bool ASMs1D::getSize (int& n1, int& n2, int& n3, int basis) const
{
n1 = this->getSize(basis);
n2 = n3 = 0;
return true;
}
int ASMs1D::getSize (int) const
{
if (!curv) return 0;

7
src/ASM/ASMs1D.h
View File

@@ -269,6 +269,13 @@ public:
//! \param[out] p3 Order in third (w) direction (always zero)
virtual bool getOrder(int& p1, int& p2, int& p3) const;
//! \brief Returns the number of nodal points in each parameter direction.
//! \param[out] n1 Number of nodes in first (u) direction
//! \param[out] n2 Number of nodes in second (v) direction
//! \param[out] n3 Number of nodes in third (w) direction
//! \param[in] basis Which basis to return size parameters for (mixed methods)
virtual bool getSize(int& n1, int& n2, int& n3, int basis) const;
protected:
Go::SplineCurve* curv; //!< Pointer to the actual spline curve object
};

7
src/ASM/ASMs2D.C
View File

@@ -1198,6 +1198,13 @@ bool ASMs2D::getOrder (int& p1, int& p2, int& p3) const
}
bool ASMs2D::getSize (int& n1, int& n2, int& n3, int basis) const
{
n3 = 0;
return this->getSize(n1,n2,basis);
}
bool ASMs2D::getSize (int& n1, int& n2, int) const
{
if (!surf) return false;

7
src/ASM/ASMs2D.h
View File

@@ -472,6 +472,13 @@ public:
//! \param[in] basis Which basis to return size parameters for (mixed methods)
virtual bool getSize(int& n1, int& n2, int basis = 0) const;
//! \brief Returns the number of nodal points in each parameter direction.
//! \param[out] n1 Number of nodes in first (u) direction
//! \param[out] n2 Number of nodes in second (v) direction
//! \param[out] n3 Number of nodes in third (w) direction
//! \param[in] basis Which basis to return size parameters for (mixed methods)
virtual bool getSize(int& n1, int& n2, int& n3, int basis) const;
//! \brief Returns the number of elements on a boundary.
virtual size_t getNoBoundaryElms(char lIndex, char ldim) const;

7
src/ASM/ASMstruct.h
View File

@@ -48,6 +48,13 @@ public:
//! \brief Resets the global element and node counters.
static void resetNumbering(int nnod = 0);
//! \brief Returns the number of nodal points in each parameter direction.
//! \param[out] n1 Number of nodes in first (u) direction
//! \param[out] n2 Number of nodes in second (v) direction
//! \param[out] n3 Number of nodes in third (w) direction
//! \param[in] basis Which basis to return size parameters for (mixed methods)
virtual bool getSize(int& n1, int& n2, int& n3, int basis = 0) const = 0;
//! \brief Projects the secondary solution field onto the primary basis.
//! \param[in] integr Object with problem-specific data and methods
virtual Go::GeomObject* evalSolution(const IntegrandBase& integr) const = 0;

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

@@ -44,17 +44,15 @@
ASMu3D::ASMu3D (unsigned char n_f)
: ASMunstruct(3,3,n_f), lrspline(0), tensorspline(0), workingEl(-1)
: ASMunstruct(3,3,n_f), lrspline(0), tensorspline(0)
{
ASMunstruct::resetNumbering(); // Replace this when going multi-patch...
swapW = false;
}
ASMu3D::ASMu3D (const ASMu3D& patch, unsigned char n_f)
: ASMunstruct(patch,n_f), lrspline(patch.lrspline), tensorspline(0), workingEl(-1)
: ASMunstruct(patch,n_f), lrspline(patch.lrspline), tensorspline(0)
{
swapW = patch.swapW;
}
size_t ASMu3D::getNodeIndex (int globalNum, bool noAddedNodes) const
@@ -63,7 +61,7 @@ size_t ASMu3D::getNodeIndex (int globalNum, bool noAddedNodes) const
if (it == MLGN.end()) return 0;
size_t inod = 1 + (it-MLGN.begin());
if (noAddedNodes && !xnMap.empty() )
if (noAddedNodes && !xnMap.empty())
{
std::map<size_t,size_t>::const_iterator it = xnMap.find(inod);
if (it != xnMap.end()) return it->second;
@@ -172,34 +170,17 @@ void ASMu3D::clear (bool retainGeometry)
nxMap.clear();
}
size_t ASMu3D::getNoNodes (int basis) const
{
return lrspline->nBasisFunctions();
}
bool ASMu3D::checkRightHandSystem ()
{
if (!lrspline || shareFE) return false;
std::cerr << "ASMu3D::checkRightHandSystem() not properly implemented yet" << std::endl;
exit(776654);
// Evaluate the spline volume at its center
double u = 0.5*(lrspline->startparam(0) + lrspline->endparam(0));
double v = 0.5*(lrspline->startparam(1) + lrspline->endparam(1));
double w = 0.5*(lrspline->startparam(2) + lrspline->endparam(2));
std::vector<Go::Point> pts;
lrspline->point(pts,u,v,w,1);
// Check that |J| = (dXdu x dXdv) * dXdw > 0.0
if ((pts[1] % pts[2]) * pts[3] > 0.0) return false;
// This patch has a negative Jacobian determinant. Probably it is modelled
// in a left-hand-system. Swap the w-parameter direction to correct for this.
// lrspline->reverseParameterDirection(2);
std::cout <<"\tSwapped."<< std::endl;
return swapW = true;
std::cout <<"ASMu3D::checkRightHandSystem(): Not available for LR-splines (ignored)." << std::endl;
return false;
}
@@ -349,18 +330,7 @@ bool ASMu3D::generateFEMTopology ()
bool ASMu3D::connectPatch (int face, ASMu3D& neighbor, int nface, int norient)
{
std::cerr << "ASMu3D::connectPatch(...) is not properly implemented yet :(" << std::endl;
exit(776654);
#if 0
if (swapW && face > 4) // Account for swapped parameter direction
face = 11-face;
if (neighbor.swapW && face > 4) // Account for swapped parameter direction
nface = 11-nface;
return this->connectBasis(face,neighbor,nface,norient);
#endif
return false;
return this->connectBasis(face,neighbor,nface,norient);
}
@@ -537,7 +507,7 @@ void ASMu3D::closeFaces (int dir, int basis, int master)
this->makePeriodic(master,master+n1*n2*(n3-1));
break;
}
#endif
#endif
}
/*!
@@ -549,9 +519,6 @@ void ASMu3D::closeFaces (int dir, int basis, int master)
void ASMu3D::constrainFace (int dir, bool open, int dof, int code)
{
if (swapW) // Account for swapped parameter direction
if (dir == 3 || dir == -3) dir = -dir;
if(open)
std::cerr << "\nWARNING: ASMu3D::constrainFace, open boundary conditions not supported yet. Treating it as closed" << std::endl;
@@ -598,9 +565,6 @@ void ASMu3D::constrainEdge (int lEdge, bool open, int dof, int code)
if(open)
std::cerr << "\nWARNING: ASMu3D::constrainEdge, open boundary conditions not supported yet. Treating it as closed" << std::endl;
if (swapW && lEdge <= 8) // Account for swapped parameter direction
lEdge += (lEdge-1)%4 < 2 ? 2 : -2;
// lEdge = 1-4, running index is u (vmin,wmin), (vmax,wmin), (vmin,wmax), (vmax,wmax)
// lEdge = 5-8, running index is v (umin,wmin), (umax,wmin), (umin,wmax), (umax,wmax)
// lEdge = 9-12, running index is w
@@ -645,15 +609,12 @@ void ASMu3D::constrainLine (int fdir, int ldir, double xi, int dof, int code)
{
std::cerr << "ASMu3D::constrainLine not implemented properly yet" << std::endl;
exit(776654);
#if 0
#if 0
if (xi < 0.0 || xi > 1.0) return;
int n1, n2, n3, node = 1;
if (!this->getSize(n1,n2,n3,1)) return;
if (swapW) // Account for swapped parameter direction
if (fdir == 3 || fdir == -3) fdir = -fdir;
switch (fdir)
{
case 1: // Right face (positive I-direction)
@@ -662,7 +623,7 @@ void ASMu3D::constrainLine (int fdir, int ldir, double xi, int dof, int code)
if (ldir == 2)
{
// Line goes in J-direction
node += n1*n2*int(0.5+(n3-1)*(swapW ? 1.0-xi : xi));
node += n1*n2*int(0.5+(n3-1)*xi);
for (int i2 = 1; i2 <= n2; i2++, node += n1)
this->prescribe(node,dof,code);
}
@@ -681,7 +642,7 @@ void ASMu3D::constrainLine (int fdir, int ldir, double xi, int dof, int code)
if (ldir == 1)
{
// Line goes in I-direction
node += n1*n2*int(0.5+(n3-1)*(swapW ? 1.0-xi : xi));
node += n1*n2*int(0.5+(n3-1)*xi);
for (int i1 = 1; i1 <= n1; i1++, node++)
this->prescribe(node,dof,code);
}
@@ -713,7 +674,7 @@ void ASMu3D::constrainLine (int fdir, int ldir, double xi, int dof, int code)
}
break;
}
#endif
#endif
}
@@ -721,20 +682,17 @@ void ASMu3D::constrainCorner (int I, int J, int K, int dof, int code)
{
std::cerr << "ASMu3D::constrainCorner not implemented properly yet" << std::endl;
exit(776654);
#if 0
#if 0
int n1, n2, n3;
if (!this->getSize(n1,n2,n3,1)) return;
if (swapW) // Account for swapped parameter direction
K = -K;
int node = 1;
if (I > 0) node += n1-1;
if (J > 0) node += n1*(n2-1);
if (K > 0) node += n1*n2*(n3-1);
this->prescribe(node,dof,code);
#endif
#endif
}
@@ -747,9 +705,6 @@ void ASMu3D::constrainNode (double xi, double eta, double zeta,
if (eta < 0.0 || eta > 1.0) return;
if (zeta < 0.0 || zeta > 1.0) return;
if (swapW) // Account for swapped parameter direction
zeta = 1.0-zeta;
#if 0
int n1, n2, n3;
if (!this->getSize(n1,n2,n3,1)) return;
@@ -777,7 +732,7 @@ bool ASMu3D::updateDirichlet (const std::map<int,RealFunc*>& func,
// 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
#if 0
std::map<int,RealFunc*>::const_iterator fit;
std::map<int,VecFunc*>::const_iterator vfit;
std::vector<DirichletFace>::const_iterator dit;
@@ -832,7 +787,7 @@ bool ASMu3D::updateDirichlet (const std::map<int,RealFunc*>& func,
// The parent class method takes care of the corner nodes with direct
// evaluation of the Dirichlet functions (since they are interpolatory)
#endif
#endif
}
#define DERR -999.99
@@ -976,24 +931,17 @@ bool ASMu3D::getOrder (int& p1, int& p2, int& p3) const
}
#if 0
bool ASMu3D::getSize (int& n1, int& n2, int& n3, int) const
{
if (!lrspline) return false;
n1 = lrspline->nBasisFunctions();
n2 = lrspline->nBasisFunctions();
n3 = lrspline->nBasisFunctions();
return true;
}
size_t ASMu3D::getNoBoundaryElms (char lIndex, char ldim) const
{
if (!lrspline) return 0;
if (ldim < 1 && lIndex > 0)
return 1;
std::cerr <<"ASMu3D::getNoBoundaryElms not implemented properly yet" << std::endl;
exit(776654);
#if 0
else if (ldim < 2 && lIndex > 0 && lIndex <= 12)
return lrspline->numCoefs((lIndex-1)/4) - lrspline->order((lIndex-1)/4) + 1;
@@ -1013,10 +961,10 @@ size_t ASMu3D::getNoBoundaryElms (char lIndex, char ldim) const
case 6:
return n1*n2;
}
#endif
return 0;
}
#endif
void ASMu3D::getGaussPointParameters (RealArray& uGP, int dir, int nGauss,
@@ -1315,7 +1263,7 @@ bool ASMu3D::integrate (Integrand& integrand,
{
std::cerr << "Haven't really figured out what this part does yet\n";
exit(42142);
#if 0
#if 0
// --- Selective reduced integration loop ----------------------------
int ip = (((i3-p3)*nRed*nel2 + i2-p2)*nRed*nel1 + i1-p1)*nRed;
@@ -1347,7 +1295,7 @@ bool ASMu3D::integrate (Integrand& integrand,
if (!integrand.reducedInt(*A,fe,X))
ok = false;
}
#endif
#endif
}
@@ -1990,8 +1938,8 @@ bool ASMu3D::evalSolution (Matrix& sField, const Vector& locSol,
for (size_t i = 0; i < nPoints; i++)
{
// fetch element containing evaluation point
// int iel = i/nPtsPerElement; // points are always listed in the same order as the elemnts
int iel = (workingEl>=0) ? workingEl : lrspline->getElementContaining(gpar[0][i], gpar[1][i], gpar[2][i]); // sadly, they are not always ordered in the same way as the elements
// sadly, points are not always ordered in the same way as the elements
int iel = lrspline->getElementContaining(gpar[0][i], gpar[1][i], gpar[2][i]);
if(iel < 0) {
std::cerr << "Logical error in evaluating results. Element at point (" << gpar[0][i] << ", "
<< gpar[1][i] << ", " << gpar[2][i] << ") not found" << std::endl;
@@ -2104,7 +2052,8 @@ bool ASMu3D::evalSolution (Matrix& sField, const IntegrandBase& integrand,
for (size_t i = 0; i < nPoints; i++)
{
// Fetch element containing evaluation point. Points are always listed
int iel = (workingEl>=0) ? workingEl : lrspline->getElementContaining(gpar[0][i], gpar[1][i], gpar[2][i]); // sadly, they are not always ordered in the same way as the elements
// sadly, points not always ordered in the same way as the elements
int iel = lrspline->getElementContaining(gpar[0][i], gpar[1][i], gpar[2][i]);
if(iel < 0) {
std::cerr << "Logical error in evaluating results. Element at point (" << gpar[0][i] << ", "

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

@@ -105,14 +105,12 @@ public:
//! \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.
//! \details If it isn't, the w-parameter direction is swapped.
bool checkRightHandSystem();
virtual bool checkRightHandSystem();
//! \brief Refines the parametrization by inserting extra knots.
//! \param[in] dir Parameter direction to refine
@@ -429,18 +427,9 @@ public:
//! \param[out] p3 Order in third (w) direction
virtual bool getOrder(int& p1, int& p2, int& p3) const;
#if 0
//! \brief Returns the number of elements on a boundary.
virtual size_t getNoBoundaryElms(char lIndex, char ldim) const;
//! \brief Returns the number of nodal points in each parameter direction.
//! \param[out] n1 Number of nodes in first (u) direction
//! \param[out] n2 Number of nodes in second (v) direction
//! \param[out] n3 Number of nodes in third (w) direction
//! \param[in] basis Which basis to return size parameters for (mixed methods)
virtual bool getSize(int& n1, int& n2, int& n3, int basis = 0) const;
#endif
// int getNodeID(size_t inod, bool = false) const { return inod; };
// unsigned char getNodalDOFs(size_t inod) const { return 0; };
// void getNoIntPoints(size_t& nPt) { };
@@ -452,19 +441,13 @@ public:
// bool injectNodeVec(const Vector& nodeVec, Vector& globVec,
// unsigned char nndof = 0) const { return false;};
private:
//! \brief Returns an index into the internal coefficient array for a node.
//! \param[in] inod 0-based node index local to current patch
int coeffInd(size_t inod) const;
protected:
LR::LRSplineVolume* lrspline; //!< Pointer to the actual spline volume object
bool swapW; //!< Has the w-parameter direction been swapped?
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; //!< The Bezier exctraction matrices for all elements
private:
Go::SplineVolume* tensorspline; //!< Pointer to original tensor spline object
@@ -472,8 +455,6 @@ private:
// 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.
mutable int workingEl; //!< here to keep track of element across function calls (to avoid topological searches all the time)
};
#endif

398
src/ASM/SAMpatchPara.C
View File

@@ -14,14 +14,7 @@
#include "SAMpatchPara.h"
#include "SystemMatrix.h"
#include "LinAlgInit.h"
#include "ASMbase.h"
#include "ASMs1D.h"
#include "ASMs2D.h"
#include "ASMs3D.h"
#include "MPC.h"
#include "GoTools/geometry/SplineCurve.h"
#include "GoTools/geometry/SplineSurface.h"
#include "GoTools/trivariate/SplineVolume.h"
#include "ASMstruct.h"
#ifdef HAS_PETSC
#include "PETScBlockMatrix.h"
@@ -69,28 +62,6 @@ bool SAMpatchPara::init (const ASMVec& model, int numNod)
}
bool SAMpatchPara::getEqns(IntVec& eqns, int f1, int f2) const
{
int nf = ndof/nnod;
int nbf = f2-f1+1;
int nbdof = nnod*nbf;
if (f1 > f2 || f1 < 0 || f2 > nf-1)
return false;
eqns.resize(nbdof);
for (int n = 0;n < nnod;n++) {
int dof = nnod*nf + f1;
int bdof = nnod*nbf;
for (int k = 0;k < nbf;k++, dof++, bdof++)
eqns[bdof] = meqn[dof];
}
return true;
}
bool SAMpatchPara::getNoDofCouplings (int ifirst, int ilast,
IntVec& d_nnz, IntVec& o_nnz) const
{
@@ -578,125 +549,32 @@ Real SAMpatchPara::normInf (const Vector& x, size_t& comp, char dofType) const
}
bool SAMpatchPara::initConstraintEqs (const std::vector<ASMbase*>& model)
bool SAMpatchPara::initConstraintEqs (const ASMVec& model)
{
// TODO: Rewrite this calling the parent-class method, and then replacing
// the mpmceq array into the ndof-sized array used here (why is it needed?)
if (!this->SAMpatch::initConstraintEqs(model))
return false;
// Estimate the size of the constraint equation array
size_t j;
MPCIter cit;
for (j = 0; j < model.size(); j++)
for (cit = model[j]->begin_MPC(); cit != model[j]->end_MPC(); cit++)
nmmceq += 1 + (*cit)->getNoMaster();
// Replace the mpmceq array by an equivalent ndof-sized array (why?)
int* new_mpmceq = new int[ndof];
memset(new_mpmceq,0,ndof*sizeof(int));
// Initialize the constraint equation arrays
mpmceq = new int[ndof];
memset(mpmceq,0,ndof*sizeof(int));
if (nceq < 1) return true;
mmceq = new int[nmmceq];
ttcc = new Real[nmmceq];
int ip = 1;
for (j = 0; j < model.size(); j++)
for (cit = model[j]->begin_MPC(); cit != model[j]->end_MPC(); cit++, ip++)
{
// Slave dof ...
int idof = madof[(*cit)->getSlave().node-1] + (*cit)->getSlave().dof - 1;
if (msc[idof-1] == 0)
if (nceq > 0)
{
int ip = 0;
MPCIter cit;
for (size_t j = 0; j < model.size(); j++)
for (cit = model[j]->begin_MPC(); cit != model[j]->end_MPC(); ++cit, ip++)
{
std::cerr <<"SAM: Ignoring constraint equation for dof "
<< idof <<" ("<< (*cit)->getSlave()
<<").\n This dof is already marked as FIXED."<< std::endl;
ip--;
nceq--;
continue;
}
else if (msc[idof-1] < 0)
{
std::cerr <<"SAM: Ignoring constraint equation for dof "
<< idof <<" ("<< (*cit)->getSlave()
<<").\n This dof is already marked as SLAVE."<< std::endl;
ip--;
nceq--;
continue;
// Slave dof ...
int idof = madof[(*cit)->getSlave().node-1] + (*cit)->getSlave().dof-2;
(*cit)->iceq = idof; // index into mpmceq for this MPC equation
new_mpmceq[idof] = mpmceq[ip];
}
}
mpmceq[idof-1] = ip;
int ipslv = ip - 1;
mmceq[ipslv] = idof;
ttcc[ipslv] = (*cit)->getSlave().coeff;
msc[idof-1] = -ip;
// Master dofs ...
for (size_t i = 0; i < (*cit)->getNoMaster(); i++)
{
idof = madof[(*cit)->getMaster(i).node-1] + (*cit)->getMaster(i).dof-1;
if (msc[idof-1] > 0)
{
int ipmst = (mpmceq[idof]++) - 1;
mmceq[ipmst] = idof;
ttcc[ipmst] = (*cit)->getMaster(i).coeff;
}
else if (msc[idof-1] < 0)
{
// Master dof is constrained (unresolved chaining)
std::cerr <<"SAM: Chained MPCs detected"
<<", slave "<< (*cit)->getSlave()
<<", master "<< (*cit)->getMaster(i)
<<" (ignored)."<< std::endl;
mpmceq[idof] = mpmceq[idof-1];
ip--;
nceq--;
break;
}
}
}
// Reset the negative values in msc before calling SYSEQ
for (ip = 0; ip < ndof; ip++)
if (msc[ip] < 0) msc[ip] = 0;
return true;
}
bool SAMpatchPara::updateConstraintEqs (const std::vector<ASMbase*>& model,
const Vector* prevSol)
{
if (nceq < 1) return true; // No constraints in this model
MPCIter cit;
for (size_t j = 0; j < model.size(); j++)
for (cit = model[j]->begin_MPC(); cit != model[j]->end_MPC(); cit++)
{
// Slave dof ...
int idof = madof[(*cit)->getSlave().node-1] + (*cit)->getSlave().dof - 1;
int ipeq = mpmceq[idof-1] - 1;
if (msc[idof-1] > 0 || mmceq[ipeq] != idof)
{
std::cerr <<" *** Corrupted SAM arrays detected in update."<< std::endl;
return false;
}
else if (!prevSol)
ttcc[ipeq] = 0.0;
else if (idof <= (int)prevSol->size())
ttcc[ipeq] = (*cit)->getSlave().coeff - (*prevSol)(idof);
else
ttcc[ipeq] = (*cit)->getSlave().coeff;
// Master dofs ...
for (size_t i = 0; prevSol && i < (*cit)->getNoMaster(); i++)
{
idof = madof[(*cit)->getMaster(i).node-1] + (*cit)->getMaster(i).dof-1;
if (msc[idof-1] > 0 && mmceq[++ipeq] == idof)
ttcc[ipeq] = (*cit)->getMaster(i).coeff;
}
}
// Swap the mpmceq array with the new one
delete[] mpmceq;
mpmceq = new_mpmceq;
return true;
}
@@ -828,7 +706,7 @@ bool SAMpatchPara::getLocalSubdomains (PetscIntMat& locSubds,
MPI_Recv(&minNodeId[0],1,MPI_INT,myRank-1,101,PETSC_COMM_WORLD,&status);
minNodeId[0]++;
}
#endif
#endif
switch (nsd) {
case 1:
@@ -902,7 +780,7 @@ bool SAMpatchPara::getLocalSubdomainsBlock (PetscIntMat& locSubds,
MPI_Recv(&minNodeId[0],1,MPI_INT,myRank-1,101,PETSC_COMM_WORLD,&status);
minNodeId[0]++;
}
#endif
#endif
switch (nsd) {
case 1:
@@ -999,20 +877,15 @@ bool SAMpatchPara::getLocalSubdomains1D (const IntVec& nxvec,
const IntVec& maxNodeId,
PetscIntMat& locSubds) const
{
// Define some parameters
const size_t npatch = patch.size();
if (nxvec.size() != npatch)
if (nxvec.size() != patch.size())
return false;
// Split the patches into smaller subdomains
for (size_t n = 0;n < npatch;n++) {
int nnod;
ASMs1D* pch1 = dynamic_cast<ASMs1D*>(patch[n]);
if (pch1)
nnod = pch1->getCurve()->numCoefs();
else
return false;
for (size_t n = 0; n < patch.size(); n++) {
int nnod, n2, n3;
if (!static_cast<ASMstruct*>(patch[n])->getSize(nnod,n2,n3) || n2 || n3)
return false;
int nx = nxvec[n];
int n1 = nnod/nx;
@@ -1057,21 +930,14 @@ bool SAMpatchPara::getLocalSubdomains2D (const IntVec& nxvec,
const IntVec& maxNodeId,
PetscIntMat& locSubds) const
{
// Define some parameters
const size_t npatch = patch.size();
if (nxvec.size() != npatch || nyvec.size() != npatch)
if (nxvec.size() != patch.size() || nyvec.size() != patch.size())
return false;
// Split the patches into smaller subdomains
for (size_t n = 0;n < npatch;n++) {
int nnod1, nnod2;
ASMs2D* pch2 = dynamic_cast<ASMs2D*>(patch[n]);
if (pch2) {
nnod1 = pch2->getSurface()->numCoefs_u();
nnod2 = pch2->getSurface()->numCoefs_v();
}
else
for (size_t n = 0; n < patch.size(); n++) {
int nnod1, nnod2, n3;
if (!static_cast<ASMstruct*>(patch[n])->getSize(nnod1,nnod2,n3) || n3)
return false;
int nx = nxvec[n];
@@ -1128,22 +994,16 @@ bool SAMpatchPara::getLocalSubdomains3D (const IntVec& nxvec,
const IntVec& maxNodeId,
PetscIntMat& locSubds) const
{
// Define some parameters
const size_t npatch = patch.size();
if (nxvec.size() != npatch || nyvec.size() != npatch || nzvec.size() != npatch)
if (nxvec.size() != patch.size() ||
nyvec.size() != patch.size() ||
nzvec.size() != patch.size())
return false;
// Split the patches into smaller subdomains
for (size_t n = 0;n < npatch;n++) {
for (size_t n = 0; n < patch.size(); n++) {
int nnod1, nnod2, nnod3;
ASMs3D* pch3 = dynamic_cast<ASMs3D*>(patch[n]);
if (pch3) {
nnod1 = pch3->getVolume()->numCoefs(0);
nnod2 = pch3->getVolume()->numCoefs(1);
nnod3 = pch3->getVolume()->numCoefs(2);
}
else
if (!static_cast<ASMstruct*>(patch[n])->getSize(nnod1,nnod2,nnod3))
return false;
int nx = nxvec[n];
@@ -1205,10 +1065,7 @@ bool SAMpatchPara::getLocalSubdomains3D (const IntVec& nxvec,
bool SAMpatchPara::getSubdomains1D (const IntVec& nxvec, int overlap,
PetscIntMat& subds) const
{
// Define some parameters
const size_t npatch = patch.size();
if (nxvec.size() != npatch)
if (nxvec.size() != patch.size())
return false;
// Overlap
@@ -1216,13 +1073,11 @@ bool SAMpatchPara::getSubdomains1D (const IntVec& nxvec, int overlap,
int ohigh = overlap/2 + overlap%2;
// Split the patches into smaller subdomains
for (size_t n = 0;n < npatch;n++) {
int nnod;
ASMs1D* pch1 = dynamic_cast<ASMs1D*>(patch[n]);
if (pch1)
nnod = pch1->getCurve()->numCoefs();
else
return false;
for (size_t n = 0; n < patch.size(); n++) {
int nnod, n2, n3;
if (!static_cast<ASMstruct*>(patch[n])->getSize(nnod,n2,n3) || n2 || n3)
return false;
int nx = nxvec[n];
int n1 = nnod/nx;
@@ -1263,10 +1118,7 @@ bool SAMpatchPara::getSubdomains2D (const IntVec& nxvec,
const IntVec& nyvec, int overlap,
PetscIntMat& subds) const
{
// Define some parameters
const size_t npatch = patch.size();
if (nxvec.size() != npatch || nyvec.size() != npatch)
if (nxvec.size() != patch.size() || nyvec.size() != patch.size())
return false;
// Overlap
@@ -1274,21 +1126,17 @@ bool SAMpatchPara::getSubdomains2D (const IntVec& nxvec,
int ohigh = overlap/2 + overlap%2;
// Split the patches into smaller subdomains
for (size_t n = 0;n < npatch;n++) {
int nnod1, nnod2;
ASMs2D* pch2 = dynamic_cast<ASMs2D*>(patch[n]);
if (pch2) {
nnod1 = pch2->getSurface()->numCoefs_u();
nnod2 = pch2->getSurface()->numCoefs_v();
}
else
for (size_t n = 0; n < patch.size(); n++) {
int nnod1, nnod2, n3;
if (!static_cast<ASMstruct*>(patch[n])->getSize(nnod1,nnod2,n3) || n3)
return false;
int nx = nxvec[n];
int ny = nyvec[n];
int n1 = nnod1/nx;
int n2 = nnod2/ny;
const IntVec& MLGN = patch[n]->getGlobalNodeNums();
for (int p2 = 0;p2 < ny;p2++) {
int jmin = 0;
@@ -1335,10 +1183,9 @@ bool SAMpatchPara::getSubdomains3D (const IntVec& nxvec,
const IntVec& nzvec, int overlap,
PetscIntMat& subds) const
{
// Define some parameters
const size_t npatch = patch.size();
if (nxvec.size() != npatch || nyvec.size() != npatch || nzvec.size() != npatch)
if (nxvec.size() != patch.size() ||
nyvec.size() != patch.size() ||
nzvec.size() != patch.size())
return false;
// Overlap
@@ -1346,24 +1193,19 @@ bool SAMpatchPara::getSubdomains3D (const IntVec& nxvec,
int ohigh = overlap/2 + overlap%2;
// Split the patches into smaller subdomains
for (size_t n = 0;n < npatch;n++) {
for (size_t n = 0; n < patch.size(); n++) {
int nnod1, nnod2, nnod3;
ASMs3D* pch3 = dynamic_cast<ASMs3D*>(patch[n]);
if (pch3) {
nnod1 = pch3->getVolume()->numCoefs(0);
nnod2 = pch3->getVolume()->numCoefs(1);
nnod3 = pch3->getVolume()->numCoefs(2);
}
else
if (!static_cast<ASMstruct*>(patch[n])->getSize(nnod1,nnod2,nnod3))
return false;
int nx = nxvec[n];
int ny = nyvec[n];
int nz = nzvec[n];
int n1 = nnod1/nx;
int n2 = nnod2/ny;
int n3 = nnod3/nz;
const IntVec& MLGN = patch[n]->getGlobalNodeNums();
for (int p3 = 0;p3 < nz;p3++) {
int kmin = 0;
@@ -1419,20 +1261,15 @@ bool SAMpatchPara::getLocalSubdomains1D (PetscIntMat& locSubds,
const IntVec& maxNodeId,
int f1, int f2) const
{
// Define some parameters
const size_t npatch = patch.size();
if (nxvec.size() != npatch)
if (nxvec.size() != patch.size())
return false;
// Split the patches into smaller subdomains
for (size_t n = 0;n < npatch;n++) {
int nnod;
ASMs1D* pch1 = dynamic_cast<ASMs1D*>(patch[n]);
if (pch1)
nnod = pch1->getCurve()->numCoefs();
else
return false;
for (size_t n = 0; n < patch.size(); n++) {
int nnod, n2, n3;
if (!static_cast<ASMstruct*>(patch[n])->getSize(nnod,n2,n3) || n2 || n3)
return false;
int nx = nxvec[n];
int n1 = nnod/nx;
@@ -1482,21 +1319,14 @@ bool SAMpatchPara::getLocalSubdomains2D (PetscIntMat& locSubds,
const IntVec& maxNodeId,
int f1, int f2) const
{
// Define some parameters
const size_t npatch = patch.size();
if (nxvec.size() != npatch || nyvec.size() != npatch)
if (nxvec.size() != patch.size() || nyvec.size() != patch.size())
return false;
// Split the patches into smaller subdomains
for (size_t n = 0;n < npatch;n++) {
int nnod1, nnod2;
ASMs2D* pch2 = dynamic_cast<ASMs2D*>(patch[n]);
if (pch2) {
nnod1 = pch2->getSurface()->numCoefs_u();
nnod2 = pch2->getSurface()->numCoefs_v();
}
else
for (size_t n = 0; n < patch.size(); n++) {
int nnod1, nnod2, n3;
if (!static_cast<ASMstruct*>(patch[n])->getSize(nnod1,nnod2,n3) || n3)
return false;
int nx = nxvec[n];
@@ -1559,31 +1389,25 @@ bool SAMpatchPara::getLocalSubdomains3D (PetscIntMat& locSubds,
const IntVec& maxNodeId,
int f1, int f2) const
{
// Define some parameters
const size_t npatch = patch.size();
if (nxvec.size() != npatch || nyvec.size() != npatch || nzvec.size() != npatch)
if (nxvec.size() != patch.size() ||
nyvec.size() != patch.size() ||
nzvec.size() != patch.size())
return false;
// Split the patches into smaller subdomains
for (size_t n = 0;n < npatch;n++) {
for (size_t n = 0; n < patch.size(); n++) {
int nnod1, nnod2, nnod3;
ASMs3D* pch3 = dynamic_cast<ASMs3D*>(patch[n]);
if (pch3) {
nnod1 = pch3->getVolume()->numCoefs(0);
nnod2 = pch3->getVolume()->numCoefs(1);
nnod3 = pch3->getVolume()->numCoefs(2);
}
else
if (!static_cast<ASMstruct*>(patch[n])->getSize(nnod1,nnod2,nnod3))
return false;
int nx = nxvec[n];
int ny = nyvec[n];
int nz = nzvec[n];
int n1 = nnod1/nx;
int n2 = nnod2/ny;
int n3 = nnod3/nz;
int nfield = f2-f1+1;
const IntVec& MLGN = patch[n]->getGlobalNodeNums();
@@ -1641,24 +1465,19 @@ bool SAMpatchPara::getSubdomains1D (PetscIntMat& subds,
const IntVec& nxvec,
int overlap, int f1, int f2) const
{
// Define some parameters
const size_t npatch = patch.size();
if (nxvec.size() != npatch)
if (nxvec.size() != patch.size())
return false;
// Overlap
int olow = overlap/2;
int ohigh = overlap/2 + overlap%2;
// Split the patches into smaller subdomains
for (size_t n = 0;n < npatch;n++) {
int nnod;
ASMs1D* pch1 = dynamic_cast<ASMs1D*>(patch[n]);
if (pch1)
nnod = pch1->getCurve()->numCoefs();
else
return false;
for (size_t n = 0; n < patch.size(); n++) {
int nnod, n2, n3;
if (!static_cast<ASMstruct*>(patch[n])->getSize(nnod,n2,n3) || n2 || n3)
return false;
int nx = nxvec[n];
int n1 = nnod/nx;
@@ -1703,10 +1522,7 @@ bool SAMpatchPara::getSubdomains2D (PetscIntMat& subds,
const IntVec& nyvec,
int overlap, int f1, int f2) const
{
// Define some parameters
const size_t npatch = patch.size();
if (nxvec.size() != npatch || nyvec.size() != npatch)
if (nxvec.size() != patch.size() || nyvec.size() != patch.size())
return false;
// Overlap
@@ -1714,21 +1530,17 @@ bool SAMpatchPara::getSubdomains2D (PetscIntMat& subds,
int ohigh = overlap/2 + overlap%2;
// Split the patches into smaller subdomains
for (size_t n = 0;n < npatch;n++) {
int nnod1, nnod2;
ASMs2D* pch2 = dynamic_cast<ASMs2D*>(patch[n]);
if (pch2) {
nnod1 = pch2->getSurface()->numCoefs_u();
nnod2 = pch2->getSurface()->numCoefs_v();
}
else
for (size_t n = 0; n < patch.size(); n++) {
int nnod1, nnod2, n3;
if (!static_cast<ASMstruct*>(patch[n])->getSize(nnod1,nnod2,n3) || n3)
return false;
int nx = nxvec[n];
int ny = nyvec[n];
int n1 = nnod1/nx;
int n2 = nnod2/ny;
int nfield = f2-f1+1;
const IntVec& MLGN = patch[n]->getGlobalNodeNums();
@@ -1783,10 +1595,9 @@ bool SAMpatchPara::getSubdomains3D (PetscIntMat& subds,
const IntVec& nzvec,
int overlap, int f1, int f2) const
{
// Define some parameters
const size_t npatch = patch.size();
if (nxvec.size() != npatch || nyvec.size() != npatch || nzvec.size() != npatch)
if (nxvec.size() != patch.size() ||
nyvec.size() != patch.size() ||
nzvec.size() != patch.size())
return false;
// Overlap
@@ -1794,24 +1605,19 @@ bool SAMpatchPara::getSubdomains3D (PetscIntMat& subds,
int ohigh = overlap/2 + overlap%2;
// Split the patches into smaller subdomains
for (size_t n = 0;n < npatch;n++) {
for (size_t n = 0; n < patch.size(); n++) {
int nnod1, nnod2, nnod3;
ASMs3D* pch3 = dynamic_cast<ASMs3D*>(patch[n]);
if (pch3) {
nnod1 = pch3->getVolume()->numCoefs(0);
nnod2 = pch3->getVolume()->numCoefs(1);
nnod3 = pch3->getVolume()->numCoefs(2);
}
else
if (!static_cast<ASMstruct*>(patch[n])->getSize(nnod1,nnod2,nnod3))
return false;
int nx = nxvec[n];
int ny = nyvec[n];
int nz = nzvec[n];
int n1 = nnod1/nx;
int n2 = nnod2/ny;
int n3 = nnod3/nz;
int nfield = f2-f1+1;
const IntVec& MLGN = patch[n]->getGlobalNodeNums();

8
src/ASM/SAMpatchPara.h
View File

@@ -53,8 +53,6 @@ public:
virtual int getMinEqNumber() const { return ieqmin; }
//! \brief Returns max global equation number for this process.
virtual int getMaxEqNumber() const { return ieqmax; }
//! \brief Returns equation numbers.
virtual bool getEqns(IntVec& eqns, int f1, int f2) const;
//! \brief Computes number of couplings for each local dof
//! in a distributed matrix.
@@ -126,12 +124,6 @@ public:
virtual bool getElmEqns(IntVec& meen, int iel, int f1, int f2,
bool globalEq = false, int nedof = 0) const;
//! \brief Updates the multi-point constraint array \a TTCC.
//! \param[in] model All spline patches in the model
//! \param[in] prevSol Previous primary solution vector in DOF-order
virtual bool updateConstraintEqs(const std::vector<ASMbase*>& model,
const Vector* prevSol = NULL);
//! \brief Expands a solution vector from equation-ordering to DOF-ordering.
//! \param[in] solVec Solution vector, length = NEQ
//! \param[out] displ Displacement vector, length = NDOF = 3*NNOD