From b169a7efecdb57eae6099fe2d1c2e2157f87053c Mon Sep 17 00:00:00 2001 From: kmo Date: Thu, 19 May 2011 11:01:15 +0000 Subject: [PATCH] Added option for element-specific debug output git-svn-id: http://svn.sintef.no/trondheim/IFEM/trunk@989 e10b68d5-8a6e-419e-a041-bce267b0401d --- .../NonlinearElasticityULMixed.C | 26 +++++++++++++++---- Apps/FiniteDefElasticity/main_NonLinEl.C | 16 +++++------- src/ASM/ASMs2Dmx.C | 14 +++++----- src/ASM/ASMs2DmxLag.C | 12 ++++++--- src/ASM/ASMs3Dmx.C | 14 +++++----- src/ASM/ASMs3DmxLag.C | 12 ++++++--- src/ASM/FiniteElement.h | 5 ++-- src/ASM/SAMpatchPara.C | 2 +- src/Utility/Utilities.C | 15 +++++++++++ src/Utility/Utilities.h | 7 +++++ 10 files changed, 86 insertions(+), 37 deletions(-) diff --git a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C index 8d5cc3f3..d13ee872 100644 --- a/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C +++ b/Apps/FiniteDefElasticity/NonlinearElasticityULMixed.C @@ -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 mixedDbgEl; //!< List of elements for additional output + bool NonlinearElasticityULMixed::initElement (const std::vector& MNPC1, const std::vector& 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; diff --git a/Apps/FiniteDefElasticity/main_NonLinEl.C b/Apps/FiniteDefElasticity/main_NonLinEl.C index 67ec3407..1859654e 100644 --- a/Apps/FiniteDefElasticity/main_NonLinEl.C +++ b/Apps/FiniteDefElasticity/main_NonLinEl.C @@ -16,6 +16,7 @@ #include "LinAlgInit.h" #include "HDF5Writer.h" #include "XMLWriter.h" +#include "Utilities.h" #include "Profiler.h" #include "VTF.h" #include @@ -23,6 +24,8 @@ #include #include +extern std::vector 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) diff --git a/src/ASM/ASMs2Dmx.C b/src/ASM/ASMs2Dmx.C index 0396be4e..a6077224 100644 --- a/src/ASM/ASMs2Dmx.C +++ b/src/ASM/ASMs2Dmx.C @@ -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; } diff --git a/src/ASM/ASMs2DmxLag.C b/src/ASM/ASMs2DmxLag.C index ac82e6c2..7877c458 100644 --- a/src/ASM/ASMs2DmxLag.C +++ b/src/ASM/ASMs2DmxLag.C @@ -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; } diff --git a/src/ASM/ASMs3Dmx.C b/src/ASM/ASMs3Dmx.C index eef0d08d..9193e839 100644 --- a/src/ASM/ASMs3Dmx.C +++ b/src/ASM/ASMs3Dmx.C @@ -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; } diff --git a/src/ASM/ASMs3DmxLag.C b/src/ASM/ASMs3DmxLag.C index c61e497c..03b7c8e3 100644 --- a/src/ASM/ASMs3DmxLag.C +++ b/src/ASM/ASMs3DmxLag.C @@ -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; } diff --git a/src/ASM/FiniteElement.h b/src/ASM/FiniteElement.h index a1320a75..12ffaa3f 100644 --- a/src/ASM/FiniteElement.h +++ b/src/ASM/FiniteElement.h @@ -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 }; diff --git a/src/ASM/SAMpatchPara.C b/src/ASM/SAMpatchPara.C index abbc566a..47586690 100644 --- a/src/ASM/SAMpatchPara.C +++ b/src/ASM/SAMpatchPara.C @@ -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; diff --git a/src/Utility/Utilities.C b/src/Utility/Utilities.C index f6941503..7c98309c 100644 --- a/src/Utility/Utilities.C +++ b/src/Utility/Utilities.C @@ -12,9 +12,24 @@ //============================================================================== #include "Utilities.h" +#include #include +void utl::parseIntegers (std::vector& 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]; diff --git a/src/Utility/Utilities.h b/src/Utility/Utilities.h index 7a60220a..f12b3712 100644 --- a/src/Utility/Utilities.h +++ b/src/Utility/Utilities.h @@ -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 :. + void parseIntegers(std::vector& 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