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

Added: static member dbgElm used to activate debug print for one element

Browse Source
This commit is contained in:
Knut Morten Okstad 2016-11-14 22:55:36 +01:00
parent 9ef54c37cf
commit c6575b6d5d
6 changed files with 154 additions and 50 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
src/ASM/ASMbase.C
View File

@ -23,6 +23,7 @@
bool ASMbase::fixHomogeneousDirichlet = true;
int ASMbase::dbgElm = 0;
/*!

2
src/ASM/ASMbase.h
View File

@ -654,6 +654,8 @@ protected:
public:
static bool fixHomogeneousDirichlet; //!< If \e true, pre-eliminate fixed DOFs
static int dbgElm; //!< One-based element index to print debugging info for
size_t idx; //!< Index of this patch in the multi-patch model
protected:

56
src/ASM/ASMs2D.C
View File

@ -688,7 +688,11 @@ bool ASMs2D::connectBasis (int edge, ASMs2D& neighbor, int nedge, bool revers,
for (i = 0; i < n1; i++, node += i1)
{
int slave = slaveNodes[revers ? n1-i-1 : i];
if (coordCheck && !neighbor.getCoord(node).equal(this->getCoord(slave),xtol))
if (!coordCheck)
ASMbase::collapseNodes(neighbor,node,*this,slave);
else if (neighbor.getCoord(node).equal(this->getCoord(slave),xtol))
ASMbase::collapseNodes(neighbor,node,*this,slave);
else
{
std::cerr <<" *** ASMs2D::connectPatch: Non-matching nodes "
<< node <<": "<< neighbor.getCoord(node)
@ -696,8 +700,6 @@ bool ASMs2D::connectBasis (int edge, ASMs2D& neighbor, int nedge, bool revers,
<< slave <<": "<< this->getCoord(slave) << std::endl;
return false;
}
else
ASMbase::collapseNodes(neighbor,node,*this,slave);
}
return true;
@ -1017,7 +1019,7 @@ size_t ASMs2D::constrainEdgeLocal (int dir, bool open, int dof, int code,
void ASMs2D::constrainCorner (int I, int J, int dof, int code, char basis)
{
int node = this->getCorner(I, J, basis);
if (node > -1)
if (node > 0)
this->prescribe(node,dof,code);
}
@ -1535,6 +1537,11 @@ bool ASMs2D::integrate (Integrand& integrand,
fe.iel = MLGE[iel];
if (fe.iel < 1) continue; // zero-area element
#ifdef SP_DEBUG
if (dbgElm < 0 && 1+iel != -dbgElm)
continue; // Skipping all elements, except for -dbgElm
#endif
int i1 = p1 + iel % nel1;
int i2 = p2 + iel / nel1;
@ -1686,10 +1693,13 @@ bool ASMs2D::integrate (Integrand& integrand,
utl::getGmat(Jac,dXidu,fe.G);
#if SP_DEBUG > 4
std::cout <<"\niel, ip = "<< iel <<" "<< ip
<<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
if (!fe.d2NdX2.empty())
std::cout <<"d2NdX2 ="<< fe.d2NdX2;
if (iel == dbgElm || iel == -dbgElm || dbgElm == 0)
{
std::cout <<"\niel, ip = "<< iel <<" "<< ip
<<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
if (!fe.d2NdX2.empty())
std::cout <<"d2NdX2 ="<< fe.d2NdX2;
}
#endif
// Cartesian coordinates of current integration point
@ -1714,6 +1724,10 @@ bool ASMs2D::integrate (Integrand& integrand,
ok = false;
A->destruct();
#ifdef SP_DEBUG
if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm
#endif
}
}
}
@ -1788,6 +1802,11 @@ bool ASMs2D::integrate (Integrand& integrand,
fe.iel = MLGE[iel];
if (fe.iel < 1) continue; // zero-area element
#ifdef SP_DEBUG
if (dbgElm < 0 && iel != -dbgElm)
continue; // Skipping all elements, except for -dbgElm
#endif
int i1 = p1 + iel % nel1;
int i2 = p2 + iel / nel1;
@ -1864,8 +1883,9 @@ bool ASMs2D::integrate (Integrand& integrand,
utl::getGmat(Jac,dXidu,fe.G);
#if SP_DEBUG > 4
std::cout <<"\niel, jp = "<< iel <<" "<< jp
<<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
if (iel == dbgElm || iel == -dbgElm || dbgElm == 0)
std::cout <<"\niel, jp = "<< iel <<" "<< jp
<<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
#endif
// Cartesian coordinates of current integration point
@ -1890,6 +1910,10 @@ bool ASMs2D::integrate (Integrand& integrand,
ok = false;
A->destruct();
#ifdef SP_DEBUG
if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm
#endif
}
}
}
@ -2134,6 +2158,11 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex,
fe.iel = abs(MLGE[doXelms+iel-1]);
if (fe.iel < 1) continue; // zero-area element
#ifdef SP_DEBUG
if (dbgElm < 0 && iel != -dbgElm)
continue; // Skipping all elements, except for -dbgElm
#endif
// Skip elements that are not on current boundary edge
bool skipMe = false;
switch (edgeDir)
@ -2220,6 +2249,11 @@ bool ASMs2D::integrate (Integrand& integrand, int lIndex,
A->destruct();
if (!ok) return false;
#ifdef SP_DEBUG
if (dbgElm < 0 && iel == -dbgElm)
break; // Skipping all elements, except for -dbgElm
#endif
}
return true;
@ -2728,7 +2762,7 @@ bool ASMs2D::evaluate (const RealFunc* func, RealArray& vec,
}
int ASMs2D::getCorner(int I, int J, int basis) const
int ASMs2D::getCorner (int I, int J, int basis) const
{
int n1, n2, node = 1;
for (char i = 1; i <= basis; i++)

60
src/ASM/ASMs3D.C
View File

@ -756,7 +756,11 @@ bool ASMs3D::connectBasis (int face, ASMs3D& neighbor, int nface, int norient,
}
int slave = slaveNodes[k][l];
if (coordCheck && !neighbor.getCoord(node).equal(this->getCoord(slave),xtol))
if (!coordCheck)
ASMbase::collapseNodes(neighbor,node,*this,slave);
else if (neighbor.getCoord(node).equal(this->getCoord(slave),xtol))
ASMbase::collapseNodes(neighbor,node,*this,slave);
else
{
std::cerr <<" *** ASMs3D::connectPatch: Non-matching nodes "
<< node <<": "<< neighbor.getCoord(node)
@ -764,8 +768,6 @@ bool ASMs3D::connectBasis (int face, ASMs3D& neighbor, int nface, int norient,
<< slave <<": "<< this->getCoord(slave) << std::endl;
return false;
}
else
ASMbase::collapseNodes(neighbor,node,*this,slave);
}
return true;
@ -1110,7 +1112,7 @@ size_t ASMs3D::constrainFaceLocal (int dir, bool open, int dof, int code,
}
std::vector<int> ASMs3D::getEdge(int lEdge, bool open, int basis) const
std::vector<int> ASMs3D::getEdge (int lEdge, bool open, int basis) const
{
std::vector<int> result;
int n1, n2, n3, n, inc = 1;
@ -1254,10 +1256,8 @@ void ASMs3D::constrainCorner (int I, int J, int K, int dof,
int code, char basis)
{
int node = this->getCorner(I, J, K, basis);
if (node < 1)
return;
this->prescribe(node,dof,code);
if (node > 0)
this->prescribe(node,dof,code);
}
@ -1785,6 +1785,11 @@ bool ASMs3D::integrate (Integrand& integrand,
fe.iel = MLGE[iel];
if (fe.iel < 1) continue; // zero-volume element
#ifdef SP_DEBUG
if (dbgElm < 0 && 1+iel != -dbgElm)
continue; // Skipping all elements, except for -dbgElm
#endif
int i1 = p1 + iel % nel1;
int i2 = p2 + (iel / nel1) % nel2;
int i3 = p3 + iel / (nel1*nel2);
@ -1946,8 +1951,13 @@ bool ASMs3D::integrate (Integrand& integrand,
utl::getGmat(Jac,dXidu,fe.G);
#if SP_DEBUG > 4
std::cout <<"\niel, ip = "<< iel <<" "<< ip
<<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
if (iel == dbgElm || iel == -dbgElm || dbgElm == 0)
{
std::cout <<"\niel, ip = "<< iel <<" "<< ip
<<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
if (!fe.d2NdX2.empty())
std::cout <<"d2NdX2 ="<< fe.d2NdX2;
}
#endif
// Cartesian coordinates of current integration point
@ -1972,6 +1982,10 @@ bool ASMs3D::integrate (Integrand& integrand,
ok = false;
A->destruct();
#ifdef SP_DEBUG
if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm
#endif
}
}
}
@ -2050,6 +2064,11 @@ bool ASMs3D::integrate (Integrand& integrand,
fe.iel = MLGE[iel];
if (fe.iel < 1) continue; // zero-volume element
#ifdef SP_DEBUG
if (dbgElm < 0 && iel != -dbgElm)
continue; // Skipping all elements, except for -dbgElm
#endif
int i1 = p1 + iel % nel1;
int i2 = p2 + (iel / nel1) % nel2;
int i3 = p3 + iel / (nel1*nel2);
@ -2129,8 +2148,9 @@ bool ASMs3D::integrate (Integrand& integrand,
utl::getGmat(Jac,dXidu,fe.G);
#if SP_DEBUG > 4
std::cout <<"\niel, jp = "<< iel <<" "<< jp
<<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
if (iel == dbgElm || iel == -dbgElm || dbgElm == 0)
std::cout <<"\niel, jp = "<< iel <<" "<< jp
<<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
#endif
// Cartesian coordinates of current integration point
@ -2155,6 +2175,10 @@ bool ASMs3D::integrate (Integrand& integrand,
ok = false;
A->destruct();
#ifdef SP_DEBUG
if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm
#endif
}
}
}
@ -2278,6 +2302,11 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex,
fe.iel = abs(MLGE[doXelms+iel]);
if (fe.iel < 1) continue; // zero-volume element
#ifdef SP_DEBUG
if (dbgElm < 0 && iel != -dbgElm)
continue; // Skipping all elements, except for -dbgElm
#endif
int i1 = p1 + iel % nel1;
int i2 = p2 + (iel / nel1) % nel2;
int i3 = p3 + iel / (nel1*nel2);
@ -2395,6 +2424,11 @@ bool ASMs3D::integrate (Integrand& integrand, int lIndex,
ok = false;
A->destruct();
#ifdef SP_DEBUG
if (dbgElm < 0 && iel == -dbgElm)
break; // Skipping all elements, except for -dbgElm
#endif
}
}
}
@ -3190,7 +3224,7 @@ bool ASMs3D::evaluate (const RealFunc* func, RealArray& vec,
}
int ASMs3D::getCorner(int I, int J, int K, int basis) const
int ASMs3D::getCorner (int I, int J, int K, int basis) const
{
int n1, n2, n3;
int node = this->findStartNode(n1,n2,n3,basis);

10
src/ASM/ASMs3Dmx.C
View File

@ -448,15 +448,14 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
// Evaluate basis function derivatives at all integration points
std::vector<std::vector<Go::BasisDerivs>> splinex(m_basis.size());
std::vector<std::vector<Go::BasisDerivs2>> splinex2(m_basis.size());
if (use2ndDer) {
if (use2ndDer)
#pragma omp parallel for schedule(static)
for (size_t i=0;i<m_basis.size();++i)
m_basis[i]->computeBasisGrid(gpar[0],gpar[1],gpar[2],splinex2[i]);
} else {
else
#pragma omp parallel for schedule(static)
for (size_t i=0;i<m_basis.size();++i)
m_basis[i]->computeBasisGrid(gpar[0],gpar[1],gpar[2],splinex[i]);
}
std::vector<size_t> elem_sizes;
for (auto& it : m_basis)
@ -553,13 +552,12 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
fe.w = gpar[2](k+1,i3-p3+1);
// Fetch basis function derivatives at current integration point
if (use2ndDer) {
if (use2ndDer)
for (size_t b = 0; b < m_basis.size(); ++b)
SplineUtils::extractBasis(splinex2[b][ip],fe.basis(b+1),dNxdu[b], d2Nxdu2[b]);
} else {
else
for (size_t b = 0; b < m_basis.size(); ++b)
SplineUtils::extractBasis(splinex[b][ip],fe.basis(b+1),dNxdu[b]);
}
// Compute Jacobian inverse of the coordinate mapping and
// basis function derivatives w.r.t. Cartesian coordinates

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

@ -637,13 +637,13 @@ void ASMu2D::constrainEdge (int dir, bool open, int dof, int code, char basis)
// get all elements and functions on this edge
std::vector<LR::Basisfunction*> edgeFunctions;
std::vector<LR::Element*> edgeElements;
lr->getEdgeFunctions(edgeFunctions, edge);
lr->getEdgeElements( edgeElements, edge);
lr->getEdgeFunctions(edgeFunctions,edge);
lr->getEdgeElements (edgeElements ,edge);
// find the corners since these are not to be included in the L2-fitting
// of the inhomogenuous dirichlet boundaries; corners are interpolatory.
// Optimization note: loop over the "edge"-container to manually pick up the
// end nodes. LRspine::getEdgeFunctions() does a global search.
// Optimization note: loop over the "edge"-container to manually pick up
// the end nodes. LRspine::getEdgeFunctions() does a global search.
std::vector<LR::Basisfunction*> c1, c2;
switch (edge)
{
@ -676,7 +676,7 @@ void ASMu2D::constrainEdge (int dir, bool open, int dof, int code, char basis)
int bcode = abs(code);
int j = 0;
for (auto b : edgeFunctions )
for (auto b : edgeFunctions)
{
de.MLGN[j++] = b->getId();
// skip corners for open boundaries
@ -993,6 +993,11 @@ bool ASMu2D::integrate (Integrand& integrand,
std::vector<LR::Element*>::iterator el = lrspline->elementBegin();
for (int iel = 1; el != lrspline->elementEnd(); el++, iel++)
{
#ifdef SP_DEBUG
if (dbgElm < 0 && iel != -dbgElm)
continue; // Skipping all elements, except for -dbgElm
#endif
FiniteElement fe(MNPC[iel-1].size());
fe.iel = MLGE[iel-1];
@ -1091,13 +1096,18 @@ bool ASMu2D::integrate (Integrand& integrand,
lrspline->computeBasis(fe.u,fe.v,spline, iel-1);
SplineUtils::extractBasis(spline,fe.N,dNdu);
#if SP_DEBUG > 4
std::cout <<"\nBasis functions at a integration point "
<<" : (u,v) = "<< spline.param[0] <<" "<< spline.param[1]
<<" left_idx = "<< spline.left_idx[0] <<" "<< spline.left_idx[1];
for (size_t ii = 0; ii < spline.basisValues.size(); ii++)
std::cout <<'\n'<< 1+ii <<'\t' << spline.basisValues[ii] <<'\t'
<< spline.basisDerivs_u[ii] <<'\t'<< spline.basisDerivs_v[ii];
std::cout << std::endl;
if (iel == dbgElm || iel == -dbgElm || dbgElm == 0)
{
std::cout <<"\nBasis functions at a integration point "
<<" : (u,v) = "<< spline.param[0] <<" "<< spline.param[1]
<<" left_idx = "<< spline.left_idx[0]
<<" "<< spline.left_idx[1];
for (size_t ii = 0; ii < spline.basisValues.size(); ii++)
std::cout <<'\n'<< 1+ii <<'\t' << spline.basisValues[ii] <<'\t'
<< spline.basisDerivs_u[ii] <<'\t'
<< spline.basisDerivs_v[ii];
std::cout << std::endl;
}
#endif
}
@ -1111,7 +1121,8 @@ bool ASMu2D::integrate (Integrand& integrand,
return false;
#if SP_DEBUG > 4
std::cout <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
if (iel == dbgElm || iel == -dbgElm || dbgElm == 0)
std::cout <<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
#endif
// Cartesian coordinates of current integration point
@ -1134,6 +1145,10 @@ bool ASMu2D::integrate (Integrand& integrand,
return false;
A->destruct();
#ifdef SP_DEBUG
if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm
#endif
}
return true;
@ -1178,6 +1193,11 @@ bool ASMu2D::integrate (Integrand& integrand,
double dA = this->getParametricArea(iel);
if (dA < 0.0) return false; // topology error (probably logic error)
#ifdef SP_DEBUG
if (dbgElm < 0 && iel != -dbgElm)
continue; // Skipping all elements, except for -dbgElm
#endif
// Set up control point (nodal) coordinates for current element
if (!this->getElementCoordinates(Xnod,iel)) return false;
@ -1230,8 +1250,9 @@ bool ASMu2D::integrate (Integrand& integrand,
return false;
#if SP_DEBUG > 4
std::cout <<"\niel, ip = "<< iel <<" "<< ip
<<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
if (iel == dbgElm || iel == -dbgElm || dbgElm == 0)
std::cout <<"\niel, ip = "<< iel <<" "<< ip
<<"\nN ="<< fe.N <<"dNdX ="<< fe.dNdX;
#endif
// Cartesian coordinates of current integration point
@ -1254,6 +1275,10 @@ bool ASMu2D::integrate (Integrand& integrand,
return false;
A->destruct();
#ifdef SP_DEBUG
if (iel == -dbgElm) break; // Skipping all elements, except for -dbgElm
#endif
}
return true;
@ -1306,6 +1331,11 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex,
std::vector<LR::Element*>::iterator el = lrspline->elementBegin();
for (int iel = 1; el != lrspline->elementEnd(); el++, iel++)
{
#ifdef SP_DEBUG
if (dbgElm < 0 && iel != -dbgElm)
continue; // Skipping all elements, except for -dbgElm
#endif
// Skip elements that are not on current boundary edge
bool skipMe = false;
switch (edgeDir)
@ -1383,6 +1413,11 @@ bool ASMu2D::integrate (Integrand& integrand, int lIndex,
return false;
A->destruct();
#ifdef SP_DEBUG
if (dbgElm < 0 && iel == -dbgElm)
break; // Skipping all elements, except for -dbgElm
#endif
}
return true;
@ -1870,9 +1905,9 @@ bool ASMu2D::updateDirichlet (const std::map<int,RealFunc*>& func,
(*mit)->setSlaveCoeff(edgeControlmatrix[dof-1][nit->first]);
else //scalar condition
(*mit)->setSlaveCoeff(edgeControlpoints[nit->first]);
#if SP_DEBUG > 1
#if SP_DEBUG > 1
std::cout <<"Updated constraint: "<< **mit;
#endif
#endif
}
}
}