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

Added option for element-specific debug output

Browse Source 
git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@989 e10b68d5-8a6e-419e-a041-bce267b0401d
This commit is contained in:
kmo 2011-05-19 11:01:15 +00:00 committed by Knut Morten Okstad
parent bf1b95f2f7
commit b169a7efec
10 changed files with 86 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C
View File

@ -14,6 +14,7 @@
#include "NonlinearElasticityULMixed.h"
#include "MaterialBase.h"
#include "FiniteElement.h"
#include "TimeDomain.h"
#include "ElmMats.h"
#include "ElmNorm.h"
#include "Utilities.h"
@ -239,17 +240,15 @@ void NonlinearElasticityULMixed::setMode (SIM::SolutionMode mode)
}
#if INT_DEBUG
static int iP = 0;
#endif
static int iP = 0; //!< Local integration point counter for debug output
std::vector<int> mixedDbgEl; //!< List of elements for additional output
bool NonlinearElasticityULMixed::initElement (const std::vector<int>& MNPC1,
const std::vector<int>& MNPC2,
size_t n1)
{
#if INT_DEBUG
iP = 0;
#endif
const size_t nen1 = MNPC1.size();
const size_t nen2 = MNPC2.size();
@ -385,6 +384,23 @@ bool NonlinearElasticityULMixed::evalIntMx (LocalIntegral*& elmInt,
if (lHaveStrains)
{
if (prm.it == 0 &&
find(mixedDbgEl.begin(),mixedDbgEl.end(),fe.iel) != mixedDbgEl.end())
{
#if INT_DEBUG > 0
if (iP == 1)
#else
if (++iP == 1)
#endif
std::cout <<"\n ** Stresses in element "<< fe.iel << std::endl;
std::cout << iP;
const RealArray& sigma = Sig;
for (size_t i = 0; i < sigma.size(); i++)
std::cout <<" "<< sigma[i];
std::cout <<" "<< Sig.vonMises() << std::endl;
}
// Compute mixed strees
Bpres = Sig.trace()/3.0;
Mpres = Press * J/Theta;

16
Apps/FiniteDefElasticity/main_NonLinEl.C
View File

@ -16,6 +16,7 @@
#include "LinAlgInit.h"
#include "HDF5Writer.h"
#include "XMLWriter.h"
#include "Utilities.h"
#include "Profiler.h"
#include "VTF.h"
#include <fstream>
@ -23,6 +24,8 @@
#include <string.h>
#include <ctype.h>
extern std::vector<int> mixedDbgEl; //!< List of elements for additional output
/*!
\brief Main program for the isogeometric finite deformation solver.
@ -131,17 +134,12 @@ int main (int argc, char** argv)
else
--i;
}
else if (!strcmp(argv[i],"-dbgElm"))
while (i < argc-1 && isdigit(argv[i+1][0]))
utl::parseIntegers(mixedDbgEl,argv[++i]);
else if (!strcmp(argv[i],"-ignore"))
while (i < argc-1 && isdigit(argv[i+1][0]))
{
char* endp = 0;
int endVal = 0;
ignoredPatches.push_back(strtol(argv[++i],&endp,10));
if (endp && *endp == ':')
endVal = strtol(endp+1,&endp,10);
while (ignoredPatches.back() < endVal)
ignoredPatches.push_back(ignoredPatches.back()+1);
}
utl::parseIntegers(ignoredPatches,argv[++i]);
else if (!strcmp(argv[i],"-vox") && i < argc-1)
VTF::vecOffset[0] = atof(argv[++i]);
else if (!strcmp(argv[i],"-voy") && i < argc-1)

14
src/ASM/ASMs2Dmx.C
View File

@ -432,7 +432,8 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
for (int i2 = p2; i2 <= n2; i2++)
for (int i1 = p1; i1 <= n1; i1++, iel++)
{
if (MLGE[iel-1] < 1) continue; // zero-area element
fe.iel = MLGE[iel-1];
if (fe.iel < 1) continue; // zero-area element
// Get element area in the parameter space
double dA = this->getParametricArea(iel);
@ -451,7 +452,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[MLGE[iel-1]-1];
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
// --- Integration loop over all Gauss points in each direction ----------
@ -493,7 +494,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand,
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,MLGE[iel-1]))
if (!glInt.assemble(elmInt,fe.iel))
return false;
}
@ -577,7 +578,8 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
for (int i2 = p2; i2 <= n2; i2++)
for (int i1 = p1; i1 <= n1; i1++, iel++)
{
if (MLGE[iel-1] < 1) continue; // zero-area element
fe.iel = MLGE[iel-1];
if (fe.iel < 1) continue; // zero-area element
// Skip elements that are not on current boundary edge
bool skipMe = false;
@ -607,7 +609,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[MLGE[iel-1]-1];
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
// --- Integration loop over all Gauss points along the edge -------------
@ -650,7 +652,7 @@ bool ASMs2Dmx::integrate (Integrand& integrand, int lIndex,
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,MLGE[iel-1]))
if (!glInt.assemble(elmInt,fe.iel))
return false;
}

12
src/ASM/ASMs2DmxLag.C
View File

@ -225,6 +225,8 @@ bool ASMs2DmxLag::integrate (Integrand& integrand,
for (int i2 = 0; i2 < nely; i2++)
for (int i1 = 0; i1 < nelx; i1++, iel++)
{
fe.iel = MLGE[iel-1];
// Set up control point coordinates for current element
if (!this->getElementCoordinates(Xnod,iel)) return false;
@ -238,7 +240,7 @@ bool ASMs2DmxLag::integrate (Integrand& integrand,
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[MLGE[iel-1]-1];
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
// --- Integration loop over all Gauss points in each direction ----------
@ -274,7 +276,7 @@ bool ASMs2DmxLag::integrate (Integrand& integrand,
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,MLGE[iel-1]))
if (!glInt.assemble(elmInt,fe.iel))
return false;
}
@ -324,6 +326,8 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, int lIndex,
for (int i2 = 0; i2 < nely; i2++)
for (int i1 = 0; i1 < nelx; i1++, iel++)
{
fe.iel = MLGE[iel-1];
// Skip elements that are not on current boundary edge
bool skipMe = false;
switch (edgeDir)
@ -348,7 +352,7 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, int lIndex,
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[MLGE[iel-1]-1];
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
// --- Integration loop over all Gauss points along the edge -------------
@ -385,7 +389,7 @@ bool ASMs2DmxLag::integrate (Integrand& integrand, int lIndex,
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,MLGE[iel-1]))
if (!glInt.assemble(elmInt,fe.iel))
return false;
}

14
src/ASM/ASMs3Dmx.C
View File

@ -468,7 +468,8 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
for (int i2 = p2; i2 <= n2; i2++)
for (int i1 = p1; i1 <= n1; i1++, iel++)
{
if (MLGE[iel-1] < 1) continue; // zero-volume element
fe.iel = MLGE[iel-1];
if (fe.iel < 1) continue; // zero-volume element
// Get element volume in the parameter space
double dV = this->getParametricVolume(iel);
@ -487,7 +488,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[MLGE[iel-1]-1];
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
// --- Integration loop over all Gauss points in each direction --------
@ -531,7 +532,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand,
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,MLGE[iel-1]))
if (!glInt.assemble(elmInt,fe.iel))
return false;
}
@ -623,7 +624,8 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
for (int i2 = p2; i2 <= n2; i2++)
for (int i1 = p1; i1 <= n1; i1++, iel++)
{
if (MLGE[iel-1] < 1) continue; // zero-volume element
fe.iel = MLGE[iel-1];
if (fe.iel < 1) continue; // zero-volume element
// Skip elements that are not on current boundary face
bool skipMe = false;
@ -665,7 +667,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[MLGE[iel-1]-1];
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
// --- Integration loop over all Gauss points in each direction --------
@ -717,7 +719,7 @@ bool ASMs3Dmx::integrate (Integrand& integrand, int lIndex,
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,MLGE[iel-1]))
if (!glInt.assemble(elmInt,fe.iel))
return false;
}

12
src/ASM/ASMs3DmxLag.C
View File

@ -251,6 +251,8 @@ bool ASMs3DmxLag::integrate (Integrand& integrand,
for (int i2 = 0; i2 < nely; i2++)
for (int i1 = 0; i1 < nelx; i1++, iel++)
{
fe.iel = MLGE[iel-1];
// Set up control point coordinates for current element
if (!this->getElementCoordinates(Xnod,iel)) return false;
@ -264,7 +266,7 @@ bool ASMs3DmxLag::integrate (Integrand& integrand,
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[MLGE[iel-1]-1];
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
// --- Integration loop over all Gauss points in each direction --------
@ -302,7 +304,7 @@ bool ASMs3DmxLag::integrate (Integrand& integrand,
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,MLGE[iel-1]))
if (!glInt.assemble(elmInt,fe.iel))
return false;
}
@ -357,6 +359,8 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, int lIndex,
for (int i2 = 0; i2 < nely; i2++)
for (int i1 = 0; i1 < nelx; i1++, iel++)
{
fe.iel = MLGE[iel-1];
// Skip elements that are not on current boundary face
bool skipMe = false;
switch (faceDir)
@ -383,7 +387,7 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, int lIndex,
// LocalIntegral pointers, of length at least the number of elements in
// the model (as defined by the highest number in the MLGE array).
// If the array is shorter than this, expect a segmentation fault.
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[MLGE[iel-1]-1];
LocalIntegral* elmInt = locInt.empty() ? 0 : locInt[fe.iel-1];
// --- Integration loop over all Gauss points in each direction --------
@ -422,7 +426,7 @@ bool ASMs3DmxLag::integrate (Integrand& integrand, int lIndex,
}
// Assembly of global system integral
if (!glInt.assemble(elmInt,MLGE[iel-1]))
if (!glInt.assemble(elmInt,fe.iel))
return false;
}

5
src/ASM/FiniteElement.h
View File

@ -25,8 +25,9 @@ class FiniteElement
{
public:
//! \brief Default constructor.
FiniteElement(size_t nb = 0) : h(0.0), N(nb), detJxW(1.0) {}
FiniteElement(size_t nb = 0) : iel(0), h(0.0), N(nb), detJxW(1.0) {}
int iel; //!< Element identifier
double u; //!< First parameter of current point
double v; //!< Second parameter of current point
double w; //!< Third parameter of current point
@ -34,7 +35,7 @@ public:
Vector N; //!< Basis function values
Vector Navg; //!< Volume-averaged basis function values
Matrix dNdX; //!< First derivatives (gradient) of the basis functions
Matrix3D d2NdX2; //!< Second derivates of the basis functions
Matrix3D d2NdX2; //!< Second derivatives of the basis functions
double detJxW; //!< Weighted determinant of the coordinate mapping
};

2
src/ASM/SAMpatchPara.C
View File

@ -196,7 +196,7 @@ bool SAMpatchPara::getElmEqns (IntVec& meen, int iel, int nedof) const
for (int j = madof[node-1]; j < madof[node]; j++)
meen.push_back(j);
}
if (meen.size() == nedof) return true;
if ((int)meen.size() == nedof) return true;
std::cerr <<"SAMpatchPara::getElmEqns: Invalid element matrix dimension "
<< nedof <<" (should have been "<< meen.size() <<")"<< std::endl;

15
src/Utility/Utilities.C
View File

@ -12,9 +12,24 @@
//==============================================================================
#include "Utilities.h"
#include <stdlib.h>
#include <ctype.h>
void utl::parseIntegers (std::vector<int>& values, const char* argv)
{
if (!argv) return;
char* endp = 0;
int endVal = 0;
values.push_back(strtol(argv,&endp,10));
if (endp && *endp == ':')
endVal = strtol(endp+1,&endp,10);
while (values.back() < endVal)
values.push_back(values.back()+1);
}
char* utl::readLine (std::istream& is)
{
static char buf[1024];

7
src/Utility/Utilities.h
View File

@ -22,6 +22,13 @@
namespace utl
{
//! \brief Parses a character string into an integer or an integer range.
//! \param values The integer value(s) is(are) appended to this vector
//! \param argv Character string with integer data
//!
//! \details An integer range is recognised through the syntax <i>:<j>.
void parseIntegers(std::vector<int>& values, const char* argv);
//! \brief Reads one line, ignoring comment lines and leading blanks.
//! \details The data read is kept in an internal static buffer.
//! \param is File stream to read from