From 71b08bd3b33270da508e899e8a2c7f78518b9414 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Wed, 5 Nov 2025 09:18:05 +0100 Subject: [PATCH 1/6] Added: ASMu2DLag::writeXML() to dump the mesh in a readable format --- src/ASM/ASMu2DLag.C | 32 ++++++++++++++++++++++++++++++++ src/ASM/ASMu2DLag.h | 3 +++ 2 files changed, 35 insertions(+) diff --git a/src/ASM/ASMu2DLag.C b/src/ASM/ASMu2DLag.C index 4d5509561..5b96c2c0e 100644 --- a/src/ASM/ASMu2DLag.C +++ b/src/ASM/ASMu2DLag.C @@ -21,6 +21,7 @@ #include "IFEM.h" #include #include +#include #ifdef USE_OPENMP #include @@ -587,3 +588,34 @@ bool ASMu2DLag::integrate (Integrand& integrand, return ok; } + + +bool ASMu2DLag::writeXML (const char* fname) const +{ + std::ofstream os(fname); + if (!os) return false; + + os <<"\n "; + for (const Vec3& X : coord) + os <<"\n "<< X; + os <<"\n "; + for (size_t nen = 1; nen <= 4; nen++) + { + bool haveElms = false; + for (const IntVec& mnpc : MNPC) + if (mnpc.size() == nen) + { + if (!haveElms) + os <<"\n "; + os <<"\n "; + for (int inod : mnpc) + os <<" "<< inod; + haveElms = true; + } + if (haveElms) + os <<"\n "; + } + os <<"\n\n"; + + return true; +} diff --git a/src/ASM/ASMu2DLag.h b/src/ASM/ASMu2DLag.h index 1e2f5be23..6b30ced43 100644 --- a/src/ASM/ASMu2DLag.h +++ b/src/ASM/ASMu2DLag.h @@ -136,6 +136,9 @@ class ASMu2DLag : public ASMs2DLag //! \note The number of element nodes must be set in \a grid on input. virtual bool tesselate(ElementBlock& grid, const int*) const; + //! \brief Dumps the mesh to the specified XML-file. + bool writeXML(const char* fname) const; + protected: //! \brief Generate thread groups using multi-coloring. void generateThreadGroupsMultiColored(bool silence, bool separateGroup1Noded); From 6600a7cd19b475258c243afa58c4fe2271a36394 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Sat, 8 Nov 2025 11:01:04 +0100 Subject: [PATCH 2/6] Changed: bool return value for utl::parseIntegers(). Added: utl::parseVec() with unit test. --- src/Utility/Test/TestUtilities.C | 81 +++++++++++++++++--------------- src/Utility/Utilities.C | 24 ++++++++-- src/Utility/Utilities.h | 11 +++-- 3 files changed, 70 insertions(+), 46 deletions(-) diff --git a/src/Utility/Test/TestUtilities.C b/src/Utility/Test/TestUtilities.C index d5fb5bd68..c28893b0a 100644 --- a/src/Utility/Test/TestUtilities.C +++ b/src/Utility/Test/TestUtilities.C @@ -10,38 +10,45 @@ //! //============================================================================== +#include "Catch2Support.h" #include "Utilities.h" #include "Vec3.h" #include "tinyxml2.h" -#include "Catch2Support.h" - #include -#include TEST_CASE("TestUtilities.ParseIntegers") { - const char* simple = "1"; - const char* ranged = "1:5"; - const char* multip = "1 3 5:8 10"; std::vector values1, values2, values3; - utl::parseIntegers(values1, simple); - utl::parseIntegers(values2, ranged); - utl::parseIntegers(values3, multip); + REQUIRE(utl::parseIntegers(values1,"1")); + REQUIRE(utl::parseIntegers(values2,"1:5")); + REQUIRE(utl::parseIntegers(values3,"1 3 5:8 10")); - int i; REQUIRE(values1.size() == 1); REQUIRE(values2.size() == 5); REQUIRE(values3.size() == 7); REQUIRE(values1.front() == 1); - for (i = 0; i < 5; i++) - REQUIRE(values2[i] == 1+i); - for (i = 0; i < 3; i++) + for (size_t i = 0; i < values2.size(); i++) + REQUIRE(values2[i] == 1+static_cast(i)); + for (int i = 0; i < 3; i++) REQUIRE(values3[i] == 1+2*i); - for (i = 2; i < 6; i++) + for (int i = 2; i < 6; i++) REQUIRE(values3[i] == 3+i); - REQUIRE(values3[6] == 10); + REQUIRE(values3.back() == 10); +} + + +TEST_CASE("TestUtilities.ParseVec3") +{ + Vec3 v1, v2, v3; + REQUIRE(utl::parseVec(v1,"1.2")); + REQUIRE(utl::parseVec(v2,"1 2")); + REQUIRE(utl::parseVec(v3,"1.2 3.4 5.6")); + + REQUIRE(v1.equal({1.2,0.0,0.0},0.0)); + REQUIRE(v2.equal({1.0,2.0,0.0},0.0)); + REQUIRE(v3.equal({1.2,3.4,5.6},0.0)); } @@ -50,7 +57,7 @@ TEST_CASE("TestUtilities.ParseKnots") std::string simple("xi 0 0.5 0.9 1.0"); strtok(const_cast(simple.c_str())," "); std::vector xi; - utl::parseKnots(xi); + REQUIRE(utl::parseKnots(xi)); REQUIRE(xi.size() == 4); REQUIRE_THAT(xi[0], WithinAbs(0.0, 1e-14)); REQUIRE_THAT(xi[1], WithinRel(0.5)); @@ -60,7 +67,7 @@ TEST_CASE("TestUtilities.ParseKnots") std::string graded5("xi C5 79 0.01 2.0"); strtok(const_cast(graded5.c_str())," "); xi.clear(); - utl::parseKnots(xi); + REQUIRE(utl::parseKnots(xi)); REQUIRE(xi.size() == 79); REQUIRE_THAT(xi[39], WithinRel(0.5)); REQUIRE_THAT(xi.front()+xi.back(), WithinRel(1.0)); @@ -118,29 +125,30 @@ TEST_CASE("TestUtilities.GetAttribute") tinyxml2::XMLDocument doc; doc.Parse(input); - REQUIRE(doc.RootElement() != nullptr); - const char* bTrue[4] = { "booltrue" , "boolone" , "boolon" , "boolyes" }; - const char* bFalse[4] = { "boolfalse", "boolzero", "booloff", "boolno" }; + const tinyxml2::XMLElement* root = doc.RootElement(); + REQUIRE(root != nullptr); - int i; - for (i = 0; i < 4; i++) { - const tinyxml2::XMLElement* elem = doc.RootElement()->FirstChildElement(bTrue[i]); + const tinyxml2::XMLElement* elem; + for (const char* bTrue : { "booltrue", "boolone", "boolon", "boolyes" }) + { + elem = root->FirstChildElement(bTrue); REQUIRE(elem != nullptr); bool b = false; REQUIRE(utl::getAttribute(elem, "bar", b)); REQUIRE(b); } - for (i = 0; i < 4; i++) { - const tinyxml2::XMLElement* elem = doc.RootElement()->FirstChildElement(bFalse[i]); + for (const char* bFalse : { "boolfalse", "boolzero", "booloff", "boolno" }) + { + elem = root->FirstChildElement(bFalse); REQUIRE(elem != nullptr); bool b = false; REQUIRE(utl::getAttribute(elem, "bar", b)); REQUIRE(!b); } - const tinyxml2::XMLElement* elem = doc.RootElement()->FirstChildElement("intval"); + elem = root->FirstChildElement("intval"); REQUIRE(elem != nullptr); int val1 = 0; REQUIRE(utl::getAttribute(elem, "bar", val1)); @@ -149,13 +157,13 @@ TEST_CASE("TestUtilities.GetAttribute") REQUIRE(utl::getAttribute(elem, "bar", val2)); REQUIRE(val2 == 1); - elem = doc.RootElement()->FirstChildElement("realval"); + elem = root->FirstChildElement("realval"); REQUIRE(elem != nullptr); double val3 = 0.0; REQUIRE(utl::getAttribute(elem, "bar", val3)); REQUIRE_THAT(val3, WithinRel(2.01)); - elem = doc.RootElement()->FirstChildElement("vecval1"); + elem = root->FirstChildElement("vecval1"); REQUIRE(elem != nullptr); Vec3 val4; REQUIRE(utl::getAttribute(elem, "bar", val4)); @@ -171,7 +179,7 @@ TEST_CASE("TestUtilities.GetAttribute") REQUIRE_THAT(val4.y, WithinRel(1.2)); REQUIRE_THAT(val4.z, WithinRel(1.2)); - elem = doc.RootElement()->FirstChildElement("vecval2"); + elem = root->FirstChildElement("vecval2"); REQUIRE(elem != nullptr); REQUIRE(utl::getAttribute(elem, "bar", val4)); REQUIRE_THAT(val4.x, WithinRel(1.2)); @@ -186,7 +194,7 @@ TEST_CASE("TestUtilities.GetAttribute") REQUIRE_THAT(val4.y, WithinRel(3.4)); REQUIRE_THAT(val4.z, WithinRel(3.4)); - elem = doc.RootElement()->FirstChildElement("vecval3"); + elem = root->FirstChildElement("vecval3"); REQUIRE(elem != nullptr); REQUIRE(utl::getAttribute(elem, "bar", val4)); REQUIRE_THAT(val4.x, WithinRel(1.2)); @@ -205,13 +213,8 @@ TEST_CASE("TestUtilities.GetAttribute") TEST_CASE("TestUtilities.GetDirs") { - int cmp1 = utl::getDirs(1); - int cmp2 = utl::getDirs(2); - int cmp3 = utl::getDirs(3); - int cmp4 = utl::getDirs(4); - - REQUIRE(cmp1 == 1); - REQUIRE(cmp2 == 12); - REQUIRE(cmp3 == 123); - REQUIRE(cmp4 == 1234); + REQUIRE(utl::getDirs(1) == 1); + REQUIRE(utl::getDirs(2) == 12); + REQUIRE(utl::getDirs(3) == 123); + REQUIRE(utl::getDirs(4) == 1234); } diff --git a/src/Utility/Utilities.C b/src/Utility/Utilities.C index 25baa530f..dac8ca72f 100644 --- a/src/Utility/Utilities.C +++ b/src/Utility/Utilities.C @@ -14,15 +14,16 @@ #include "Utilities.h" #include "Vec3.h" #include "tinyxml2.h" +#include #include #include -#include -void utl::parseIntegers (std::vector& values, const char* argv) +bool utl::parseIntegers (std::vector& values, const char* argv) { - if (!argv) return; + if (!argv) return false; + size_t old = values.size(); char* endp = const_cast(argv); while (endp && strlen(endp) > 0) { @@ -34,6 +35,23 @@ void utl::parseIntegers (std::vector& values, const char* argv) values.push_back(values.back()+1); } } + + return values.size() > old; +} + + +bool utl::parseVec (Vec3& value, const char* argv) +{ + if (!argv) return false; + + char* endp = const_cast(argv); + for (int i = 0; i < 3; i++) + if (endp && strlen(endp) > 0) + value[i] = strtod(endp,&endp); + else + value[i] = Real(0); + + return true; } diff --git a/src/Utility/Utilities.h b/src/Utility/Utilities.h index afda170b0..41b17745d 100644 --- a/src/Utility/Utilities.h +++ b/src/Utility/Utilities.h @@ -15,11 +15,9 @@ #define _UTILITIES_H #include "matrix.h" -#include -#include -#include #include #include +#include namespace tinyxml2 { class XMLElement; @@ -38,7 +36,12 @@ namespace utl //! \param[in] argv Character string with integer data //! //! \details An integer range is recognised through the syntax \a i:j. - void parseIntegers(std::vector& values, const char* argv); + bool parseIntegers(std::vector& values, const char* argv); + + //! \brief Parses a character string into a point vector. + //! \param[out] value The parsed point vector + //! \param[in] argv Character string with data + bool parseVec(Vec3& value, const char* argv); //! \brief Parses a (possibly graded) sequence of knot values. //! \param xi The knot value(s) is/are appended to this vector From 5ccee6bbc0d3bd8c5f7eef352d0116e80f923675 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Sun, 9 Nov 2025 07:25:16 +0100 Subject: [PATCH 3/6] Added: virtual method ASMbase::printElmInfo() --- src/ASM/ASMbase.h | 2 ++ src/SIM/SIMoutput.C | 2 ++ 2 files changed, 4 insertions(+) diff --git a/src/ASM/ASMbase.h b/src/ASM/ASMbase.h index 6aec0d545..604fdd234 100644 --- a/src/ASM/ASMbase.h +++ b/src/ASM/ASMbase.h @@ -317,6 +317,8 @@ class ASMbase void printNodes(std::ostream& os) const; //! \brief Prints out element connections of this patch to the given stream. void printElements(std::ostream& os) const; + //! \brief Prints out additional app-dependent element information. + virtual void printElmInfo(int, const IntegrandBase*) const {} //! \brief Increase all global node numbers by \a nshift. virtual void shiftGlobalNodeNums(int nshift); diff --git a/src/SIM/SIMoutput.C b/src/SIM/SIMoutput.C index a5988d162..5f4145704 100644 --- a/src/SIM/SIMoutput.C +++ b/src/SIM/SIMoutput.C @@ -509,6 +509,8 @@ void SIMoutput::preprocessResPtGroup (std::string& ptFile, ResPointVec& points) IFEM::cout <<", element centers will be used."<< std::endl; else IFEM::cout <<", nodal points will be used."<< std::endl; + if (pch && pt.iel > 0) + pch->printElmInfo(pt.iel,myProblem); } } From 87014c12624c8b45e645dc6f8c168349f49bc277 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Sun, 9 Nov 2025 07:40:59 +0100 Subject: [PATCH 4/6] Fixed: Missing argument to NextSiblingElement() calls. Added: Print out global element/node ID for result points if multi-patch or the global ID is not equal to the local index. Changed: Reduce scope of some temporaries in SIMoutput.C. --- src/SIM/SIM1D.C | 2 +- src/SIM/SIM2D.C | 2 +- src/SIM/SIM3D.C | 2 +- src/SIM/SIMinput.C | 2 +- src/SIM/SIMoutput.C | 85 +++++++++++++++++++++++++-------------------- src/SIM/TimeStep.C | 2 +- 6 files changed, 52 insertions(+), 43 deletions(-) diff --git a/src/SIM/SIM1D.C b/src/SIM/SIM1D.C index 203258204..ce699a0de 100644 --- a/src/SIM/SIM1D.C +++ b/src/SIM/SIM1D.C @@ -148,7 +148,7 @@ bool SIM1D::parseGeometryTag (const tinyxml2::XMLElement* elem) utl::getAttribute(elem,"offset",offset); const tinyxml2::XMLElement* child = elem->FirstChildElement("connection"); - for (; child; child = child->NextSiblingElement()) + for (; child; child = child->NextSiblingElement("connection")) { ASM::Interface ifc; if (utl::getAttribute(child,"master",ifc.master)) diff --git a/src/SIM/SIM2D.C b/src/SIM/SIM2D.C index d1094e48e..51d09171c 100644 --- a/src/SIM/SIM2D.C +++ b/src/SIM/SIM2D.C @@ -187,7 +187,7 @@ bool SIM2D::parseGeometryTag (const tinyxml2::XMLElement* elem) std::vector top; const tinyxml2::XMLElement* child = elem->FirstChildElement("connection"); - for (; child; child = child->NextSiblingElement()) + for (; child; child = child->NextSiblingElement("connection")) { ASM::Interface ifc; bool rever = false, periodic = false; diff --git a/src/SIM/SIM3D.C b/src/SIM/SIM3D.C index e37a6f1b5..fdce6c69c 100644 --- a/src/SIM/SIM3D.C +++ b/src/SIM/SIM3D.C @@ -178,7 +178,7 @@ bool SIM3D::parseGeometryTag (const tinyxml2::XMLElement* elem) utl::getAttribute(elem,"offset",offset); const tinyxml2::XMLElement* child = elem->FirstChildElement("connection"); - for (; child; child = child->NextSiblingElement()) + for (; child; child = child->NextSiblingElement("connection")) { ASM::Interface ifc; bool periodic = false; diff --git a/src/SIM/SIMinput.C b/src/SIM/SIMinput.C index 4fe9da008..d1e1727d3 100644 --- a/src/SIM/SIMinput.C +++ b/src/SIM/SIMinput.C @@ -451,7 +451,7 @@ bool SIMinput::parseBCTag (const tinyxml2::XMLElement* elem) else if (!strcasecmp(elem->Value(),"propertycodes")) { const tinyxml2::XMLElement* code = elem->FirstChildElement("code"); - for (; code; code = code->NextSiblingElement()) + for (; code; code = code->NextSiblingElement("code")) { int icode = 0; utl::getAttribute(code,"value",icode); diff --git a/src/SIM/SIMoutput.C b/src/SIM/SIMoutput.C index 5f4145704..4f27b7142 100644 --- a/src/SIM/SIMoutput.C +++ b/src/SIM/SIMoutput.C @@ -106,8 +106,7 @@ bool SIMoutput::parseOutputTag (const tinyxml2::XMLElement* elem) { IFEM::cout <<" Parsing <"<< elem->Value() <<">"<< std::endl; - const char* funcval = utl::getValue(elem,"function"); - if (funcval) + if (const char* funcval = utl::getValue(elem,"function"); funcval) { std::string name; utl::getAttribute(elem,"name",name); @@ -141,7 +140,7 @@ bool SIMoutput::parseOutputTag (const tinyxml2::XMLElement* elem) // Parse the result point specifications. // Can either be explicit points, lines or a grid. const tinyxml2::XMLElement* point = elem->FirstChildElement("point"); - for (int i = 1; point; i++, point = point->NextSiblingElement()) + for (int i = 1; point; i++, point = point->NextSiblingElement("point")) { ResultPoint thePoint(3); if (int patch = 0; utl::getAttribute(point,"patch",patch) && patch > 0) @@ -151,20 +150,22 @@ bool SIMoutput::parseOutputTag (const tinyxml2::XMLElement* elem) if (utl::getAttribute(point,"node",thePoint.inod) && thePoint.inod > 0) IFEM::cout <<" node = "<< thePoint.inod; else + { IFEM::cout <<" xi ="; - if (utl::getAttribute(point,"u",thePoint.u[0])) - IFEM::cout <<' '<< thePoint.u[0]; - if (utl::getAttribute(point,"v",thePoint.u[1])) - IFEM::cout <<' '<< thePoint.u[1]; - if (utl::getAttribute(point,"w",thePoint.u[2])) - IFEM::cout <<' '<< thePoint.u[2]; + if (utl::getAttribute(point,"u",thePoint.u[0])) + IFEM::cout <<' '<< thePoint.u[0]; + if (utl::getAttribute(point,"v",thePoint.u[1])) + IFEM::cout <<' '<< thePoint.u[1]; + if (utl::getAttribute(point,"w",thePoint.u[2])) + IFEM::cout <<' '<< thePoint.u[2]; + } IFEM::cout << std::endl; addPoint(thePoint); } const tinyxml2::XMLElement* line = elem->FirstChildElement("line"); - for (int j = 1; line; j++, line = line->NextSiblingElement()) + for (int j = 1; line; j++, line = line->NextSiblingElement("line")) { ResultPoint thePoint(3); if (int patch = 0; utl::getAttribute(line,"patch",patch) && patch > 0) @@ -206,26 +207,22 @@ bool SIMoutput::parseOutputTag (const tinyxml2::XMLElement* elem) } for (const tinyxml2::XMLElement* child = elem->FirstChildElement("elements"); - child; child = child->NextSiblingElement()) + child; child = child->NextSiblingElement("elements")) if (const tinyxml2::XMLNode* elist = child->FirstChild(); elist) { ResultPoint thePoint; if (int patch = 0; utl::getAttribute(child,"patch",patch) && patch > 0) thePoint.patch = patch; - IntVec elements; - utl::parseIntegers(elements,elist->Value()); - for (int iel : elements) - { - thePoint.iel = iel; - addPoint(thePoint); - } - - if (!elements.empty()) + if (IntVec elements; utl::parseIntegers(elements,elist->Value())) { IFEM::cout <<"\tElement center output: P"<< thePoint.patch <<" ("; - for (const ResultPoint& rp : myPoints.back().second) - IFEM::cout <<" "<< rp.iel; + for (int iel : elements) + { + IFEM::cout <<" "<< iel; + thePoint.iel = iel; + addPoint(thePoint); + } IFEM::cout <<" )"<< std::endl; } } @@ -236,7 +233,7 @@ bool SIMoutput::parseOutputTag (const tinyxml2::XMLElement* elem) else if (grid) idxGrid = myPoints.size(); - for (int g = 1; grid; g++, grid = grid->NextSiblingElement()) + for (int g = 1; grid; g++, grid = grid->NextSiblingElement("grid")) { ResultPoint thePoint; if (int patch = 0; utl::getAttribute(grid,"patch",patch) && patch > 0) @@ -413,7 +410,7 @@ void SIMoutput::preprocessResPtGroup (std::string& ptFile, ResPointVec& points) { if (points.empty()) return; - // Check if we have Cartesian grid input + // Check if we have a regular Cartesian grid input size_t nX = 0; IntMat pointMap; if (points.front().inod < 0 && points.back().inod < 0) @@ -431,8 +428,7 @@ void SIMoutput::preprocessResPtGroup (std::string& ptFile, ResPointVec& points) for (ResPointVec::iterator pit = points.begin(); pit != points.end();) { bool pointIsOK = false; - ASMbase* pch = this->getPatch(pit->patch,true); - if (pch && !pch->empty()) + if (ASMbase* pch = this->getPatch(pit->patch,true); pch && !pch->empty()) { if ((pointIsOK = pit->inod > 0)) // A nodal number is specified pit->X = pch->getCoord(pit->inod); @@ -441,8 +437,10 @@ void SIMoutput::preprocessResPtGroup (std::string& ptFile, ResPointVec& points) pointIsOK = fabs(pch->findPoint(pit->X,pit->u)) < 1.0e-3; else if (pit->npar > 0) // A parametric point is specified + { + pit->npar = pch->getNoParamDim(); pointIsOK = (pit->inod = pch->evalPoint(pit->u,pit->u,pit->X)) >= 0; - + } else if (pit->iel < 1) // All nodes or element centers pointIsOK = true; @@ -455,8 +453,6 @@ void SIMoutput::preprocessResPtGroup (std::string& ptFile, ResPointVec& points) if (pointIsOK) { - if (pit->npar > 0) - pit->npar = pch->getNoParamDim(); ++pit; ++iPoint; if (nX > 0) @@ -498,10 +494,24 @@ void SIMoutput::preprocessResPtGroup (std::string& ptFile, ResPointVec& points) IFEM::cout <<", node #"<< pt.inod; else if (pt.iel > 0) IFEM::cout <<", element #"<< pt.iel; - if (pt.inod > 0 && myModel.size() > 1) + ASMbase* pch = this->getPatch(pt.patch,true); + if (pch) { - ASMbase* pch = this->getPatch(pt.patch,true); - IFEM::cout <<", global #"<< pch->getNodeID(pt.inod); + int globalId = 0; + if (pt.inod > 0) + { + if (int nodeId = pch->getNodeID(pt.inod); + nodeId != pt.inod || myModel.size() > 1) + globalId = nodeId; + } + else if (pt.iel > 0) + { + if (int elmId = pch->getElmID(pt.iel); + elmId != pt.iel || myModel.size() > 1) + globalId = elmId; + } + if (globalId > 0) + IFEM::cout <<", global #"<< globalId; } if (pt.npar > 0 || pt.iel > 0) IFEM::cout <<", X = "<< pt.X << std::endl; @@ -541,7 +551,7 @@ void SIMoutput::preprocessResPtGroup (std::string& ptFile, ResPointVec& points) for (const ResultPoint& pt : points) { os << 0.0 <<" "; // dummy time - for (short int k = 0; k < pt.npar; k++) + for (short int k = 0; k < nsd; k++) os << std::setw(flWidth) << pt.X[k]; os << std::endl; } @@ -563,8 +573,7 @@ bool SIMoutput::merge (SIMbase* that, const std::map* old2new, for (ResultPoint& rp : rptp.second) { if (++iPoint == 1) IFEM::cout <<'\n'; - ASMbase* pch = that->getPatch(rp.patch,true); - if (pch && rp.inod > 0) + if (ASMbase* pch = that->getPatch(rp.patch,true); pch && rp.inod > 0) IFEM::cout <<"Result point #"<< iPoint <<": patch #"<< rp.patch <<", local node "<< rp.inod <<" global node "<< pch->getNodeID(rp.inod) << std::endl; @@ -1798,11 +1807,11 @@ bool SIMoutput::dumpMatlabGrid (std::ostream& os, const std::string& name, os <<"]="<< name << std::endl; // Write out nodes for the specified topology sets - TopologySet::const_iterator tit; for (const std::string& setname : sets) { std::set nodeSet; - if ((tit = myEntitys.find(setname)) != myEntitys.end()) + if (TopologySet::const_iterator tit = myEntitys.find(setname); + tit != myEntitys.end()) for (const TopItem& top : tit->second) if (ASMbase* pch = this->getPatch(top.patch); pch) if (top.idim+1 == pch->getNoParamDim()) @@ -2023,7 +2032,7 @@ bool SIMoutput::dumpResults (const Vector& psol, double time, } bool newLine = false; - // Convenience lambda returning proper newline character before identifier. + // Convenience lambda returning proper newline char before identifier. auto&& newLineChr = [&newLine]() { return newLine ? "\n\t\t" : ":\t"; }; if (formatted && !sol1.empty()) diff --git a/src/SIM/TimeStep.C b/src/SIM/TimeStep.C index 4df2b38b9..96f69aeda 100644 --- a/src/SIM/TimeStep.C +++ b/src/SIM/TimeStep.C @@ -118,7 +118,7 @@ bool TimeStep::parse (const tinyxml2::XMLElement* elem) if (f2 > 1.0) f2 = 1.0; const tinyxml2::XMLElement* child = elem->FirstChildElement("step"); - for (; child; child = child->NextSiblingElement()) + for (; child; child = child->NextSiblingElement("step")) { double start = 0.0, end = 0.0; utl::getAttribute(child,"start",start); From 8d15cc8f5acb5b2bb7e715509e1d69989c524196 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Fri, 14 Nov 2025 14:26:16 +0100 Subject: [PATCH 5/6] Added: Protected method evalSolPt() evaluating a nodal field at a point. The method is overriden in ASMu2DLag with the triangle element support, and will skip all 0D/1D elements when evaluating element center results, since these elements are assumed to not have any results. --- src/ASM/ASMs2DLag.C | 64 +++++++++++++++++++++++++++++++-------------- src/ASM/ASMs2DLag.h | 12 +++++++++ src/ASM/ASMs3DLag.C | 62 +++++++++++++++++++++++++++++-------------- src/ASM/ASMs3DLag.h | 13 +++++++++ src/ASM/ASMu2DLag.C | 30 +++++++++++++++++++++ src/ASM/ASMu2DLag.h | 12 +++++++++ 6 files changed, 154 insertions(+), 39 deletions(-) diff --git a/src/ASM/ASMs2DLag.C b/src/ASM/ASMs2DLag.C index 85553e948..7b97e18c6 100644 --- a/src/ASM/ASMs2DLag.C +++ b/src/ASM/ASMs2DLag.C @@ -718,7 +718,8 @@ bool ASMs2DLag::tesselate (ElementBlock& grid, const int*) const } -void ASMs2DLag::constrainEdge (int dir, bool open, int dof, int code, char basis) +void ASMs2DLag::constrainEdge (int dir, bool open, int dof, int code, + char basis) { this->ASMs2D::constrainEdge(dir, open, dof, code > 0 ? -code : code, basis); } @@ -741,39 +742,62 @@ bool ASMs2DLag::evalSolution (Matrix& sField, const Vector& locSol, size_t nCmp = locSol.size() / this->getNoProjectionNodes(); size_t ni = gpar ? gpar[0].size() : nel; size_t nj = gpar ? gpar[1].size() : 1; - size_t nen = p1*p2; sField.resize(nCmp,ni*nj); - Matrix elmSol(nCmp,nen); - RealArray N(nen), val; + RealArray N, val; double xi = 0.0, eta = 0.0; if (!gpar && !Lagrange::computeBasis(N,p1,xi,p2,eta)) return false; - size_t ip = 1; + size_t ip = 0; int iel = 0; for (size_t j = 0; j < nj; j++) - for (size_t i = 0; i < ni; i++, ip++) + for (size_t i = 0; i < ni; i++) { if (gpar) - { iel = this->findElement(gpar[0][i], gpar[1][j], &xi, &eta); - if (iel < 1 || iel > static_cast(nel)) - return false; - if (!Lagrange::computeBasis(N,p1,xi,p2,eta)) - return false; - } - else - iel++; + else // evaluate at element centers + --iel; - for (size_t a = 1; a <= nen; a++) - elmSol.fillColumn(a, locSol.ptr() + nCmp*MNPC[iel-1][a-1]); - - elmSol.multiply(N,val); - sField.fillColumn(ip,val); + if (!this->evalSolPt(iel,xi,eta,nCmp,locSol,val,N)) + return false; + else if (!val.empty()) // skip elements with no results + sField.fillColumn(++ip,val); } + sField.resize(nCmp,ip); + return true; +} + + +/*! + If \a iel is negative, it is assumed that the basis function values aready + have been evaluated in \a N, and the ξ and η parameters are not used. + If \a iel is positive, the basis functions are evaluated at point ξ,η. +*/ + +bool ASMs2DLag::evalSolPt (int iel, double xi, double eta, size_t nCmp, + const Vector& pchSol, RealArray& ptSol, + RealArray& N) const +{ + if (iel < 0 && -iel <= static_cast(nel)) // assume N is already evaluated + iel = -iel-1; + else if (iel < 1 || iel > static_cast(nel)) + return false; + else if (!Lagrange::computeBasis(N,p1,xi,p2,eta)) + return false; + else + --iel; + + ptSol.assign(nCmp,0.0); + for (size_t a = 0; a < N.size(); a++) + { + size_t ip = nCmp*MNPC[iel][a]; + for (size_t i = 0; i < nCmp; i++) + ptSol[i] += N[a]*pchSol[ip++]; + } + return true; } @@ -995,7 +1019,7 @@ int ASMs2DLag::findElement (double u, double v, double* xi, double* eta) const if (!surf) { std::cerr <<" *** ASMs2DLag::findElement: No spline geometry"<< std::endl; - return -1; + return 0; } int ku, kv; diff --git a/src/ASM/ASMs2DLag.h b/src/ASM/ASMs2DLag.h index a25229e47..154875031 100644 --- a/src/ASM/ASMs2DLag.h +++ b/src/ASM/ASMs2DLag.h @@ -290,6 +290,18 @@ class ASMs2DLag : public ASMs2D, protected ASMLagBase static void createMNPC(size_t nx, size_t ny, int p1, int p2, IntMat& MNPC); protected: + //! \brief Evaluates a nodal solution field at specified point in an element. + //! \param[in] iel 1-based element index local to current patch + //! \param[in] xi First parametric coordinate, ξ ∈ [-1,1] + //! \param[in] eta Second parametric coordinate, η ∈ [-1,1] + //! \param[in] nCmp Number of solution components + //! \param[in] pchSol Nodal values of the field to evaluate + //! \param[out] ptSol Solution field values at the specified point + //! \param[out] N Basis function values at the sepcified point + virtual bool evalSolPt(int iel, double xi, double eta, + size_t nCmp, const Vector& pchSol, + RealArray& ptSol, RealArray& N) const; + size_t nx; //!< Number of nodes in first parameter direction size_t ny; //!< Number of nodes in second parameter direction int p1; //!< Polynomial order in first parameter direction diff --git a/src/ASM/ASMs3DLag.C b/src/ASM/ASMs3DLag.C index 86671dfcb..e89a7f001 100644 --- a/src/ASM/ASMs3DLag.C +++ b/src/ASM/ASMs3DLag.C @@ -952,7 +952,7 @@ int ASMs3DLag::findElement (double u, double v, double w, if (!svol) { std::cerr <<" *** ASMs3DLag::findElement: No spline geometry"<< std::endl; - return -1; + return 0; } std::array,3> knot; @@ -1001,41 +1001,65 @@ bool ASMs3DLag::evalSolution (Matrix& sField, const Vector& locSol, size_t ni = gpar ? gpar[0].size() : nel; size_t nj = gpar ? gpar[1].size() : 1; size_t nk = gpar ? gpar[2].size() : 1; - size_t nen = p1*p2*p3; sField.resize(nCmp,ni*nj*nk); - Matrix elmSol(nCmp,nen); - RealArray N(nen), val; + RealArray N, val; double xi = 0.0, eta = 0.0, zeta = 0.0; if (!gpar && !Lagrange::computeBasis(N,p1,xi,p2,eta,p3,zeta)) return false; - size_t ip = 1; + size_t ip = 0; int iel = 0; for (size_t k = 0; k < nk; k++) for (size_t j = 0; j < nj; j++) - for (size_t i = 0; i < ni; i++, ip++) + for (size_t i = 0; i < ni; i++) { if (gpar) - { iel = this->findElement(gpar[0][i], gpar[1][j], gpar[2][k], &xi, &eta, &zeta); - if (iel < 1 || iel > static_cast(nel)) - return false; - if (!Lagrange::computeBasis(N,p1,xi,p2,eta,p3,zeta)) - return false; - } - else - iel++; + else // evaluate at element centers + --iel; - for (size_t a = 1; a <= nen; a++) - elmSol.fillColumn(a, locSol.ptr() + nCmp*MNPC[iel-1][a-1]); - - elmSol.multiply(N,val); - sField.fillColumn(ip,val); + if (!this->evalSolPt(iel,xi,eta,zeta,nCmp,locSol,val,N)) + return false; + else if (!val.empty()) // skip elements with no results + sField.fillColumn(++ip,val); } + sField.resize(nCmp,ip); + return true; +} + + +/*! + If \a iel is negative, it is assumed that the basis function values aready + have been evaluated in \a N, and the ξ, η and ζ parameters are + not used. If \a iel is positive, the basis functions are evaluated at the + point ξ,η,ζ. +*/ + +bool ASMs3DLag::evalSolPt (int iel, double xi, double eta, double zeta, + size_t nCmp, const Vector& pchSol, RealArray& ptSol, + RealArray& N) const +{ + if (iel < 0 && -iel <= static_cast(nel)) // assume N is already evaluated + iel = -iel-1; + else if (iel < 1 || iel > static_cast(nel)) + return false; + else if (!Lagrange::computeBasis(N,p1,xi,p2,eta,p3,zeta)) + return false; + else + --iel; + + ptSol.assign(nCmp,0.0); + for (size_t a = 0; a < N.size(); a++) + { + size_t ip = nCmp*MNPC[iel][a]; + for (size_t i = 0; i < nCmp; i++) + ptSol[i] += N[a]*pchSol[ip++]; + } + return true; } diff --git a/src/ASM/ASMs3DLag.h b/src/ASM/ASMs3DLag.h index c49e5695e..30b36f7a7 100644 --- a/src/ASM/ASMs3DLag.h +++ b/src/ASM/ASMs3DLag.h @@ -313,6 +313,19 @@ class ASMs3DLag : public ASMs3D, protected ASMLagBase int p1, int p2, int p3, IntMat& MNPC); protected: + //! \brief Evaluates a nodal solution field at specified point in an element. + //! \param[in] iel 1-based element index local to current patch + //! \param[in] xi First parametric coordinate in range [-1,1] + //! \param[in] eta Second parametric coordinate in range [-1,1] + //! \param[in] zeta Second parametric coordinate in range [-1,1] + //! \param[in] nCmp Number of solution components + //! \param[in] pchSol Nodal values of the field to evaluate + //! \param[out] ptSol Solution field values at the specified point + //! \param[out] N Basis function values at the sepcified point + virtual bool evalSolPt(int iel, double xi, double eta, double zeta, + size_t nCmp, const Vector& pchSol, + RealArray& ptSol, RealArray& N) const; + size_t nx; //!< Number of nodes in first parameter direction size_t ny; //!< Number of nodes in second parameter direction size_t nz; //!< Number of nodes in third parameter direction diff --git a/src/ASM/ASMu2DLag.C b/src/ASM/ASMu2DLag.C index 5b96c2c0e..81370a5ec 100644 --- a/src/ASM/ASMu2DLag.C +++ b/src/ASM/ASMu2DLag.C @@ -590,6 +590,36 @@ bool ASMu2DLag::integrate (Integrand& integrand, } +/*! + If \a iel is negative, it is assumed that the basis function values aready + have been evaluated in \a N, and the ξ and η parameters are not used. + If \a iel is positive, the basis functions are evaluated at point ξ,η. + The 4-noded elements are handled by forwarding to the parent class method. + The 3-noded elements (linear triangles) are handled separately by setting up + a temporary basis function array \a N3 and not touching the incoming \a N. + This method will ignore elements with less than 3 nodes (\a ptSol is empty). +*/ + +bool ASMu2DLag::evalSolPt (int iel, double xi, double eta, size_t nCmp, + const Vector& pchSol, RealArray& ptSol, + RealArray& N) const +{ + ptSol.clear(); + size_t jel = iel > 0 ? iel : -iel; + if (jel < 1 || jel > MNPC.size()) + return false; + else if (MNPC[--jel].size() > 3) // 4-noded 2D element + return this->ASMs2DLag::evalSolPt(iel,xi,eta,nCmp,pchSol,ptSol,N); + else if (MNPC[jel].size() < 3) + return true; // no results for 0D and 1D elements + + // This element is a linear triangle + RealArray N3(3,1.0/3.0); + if (iel > 0) N3 = { xi, eta, 1.0-xi-eta }; + return this->ASMs2DLag::evalSolPt(-1-jel,xi,eta,nCmp,pchSol,ptSol,N3); +} + + bool ASMu2DLag::writeXML (const char* fname) const { std::ofstream os(fname); diff --git a/src/ASM/ASMu2DLag.h b/src/ASM/ASMu2DLag.h index 6b30ced43..0d1dc9dc2 100644 --- a/src/ASM/ASMu2DLag.h +++ b/src/ASM/ASMu2DLag.h @@ -140,6 +140,18 @@ class ASMu2DLag : public ASMs2DLag bool writeXML(const char* fname) const; protected: + //! \brief Evaluates a nodal solution field at specified point in an element. + //! \param[in] iel 1-based element index local to current patch + //! \param[in] xi First parametric coordinate in range [-1,1] + //! \param[in] eta Second parametric coordinate in range [-1,1] + //! \param[in] nCmp Number of solution components + //! \param[in] pchSol Nodal values of the field to evaluate + //! \param[out] ptSol Solution field values at the specified point + //! \param[out] N Basis function values at the sepcified point + virtual bool evalSolPt(int iel, double xi, double eta, + size_t nCmp, const Vector& pchSol, + RealArray& ptSol, RealArray& N) const; + //! \brief Generate thread groups using multi-coloring. void generateThreadGroupsMultiColored(bool silence, bool separateGroup1Noded); From 077b58c3748ff6d0b55989b8cc507c36203a960f Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Tue, 18 Nov 2025 14:09:39 +0100 Subject: [PATCH 6/6] Added: Print out the secondary solution component names. Added: Specification of which solution components to output. --- src/SIM/SIMoutput.C | 38 ++++++++++++++++++++++++++++++++------ src/SIM/SIMoutput.h | 2 ++ 2 files changed, 34 insertions(+), 6 deletions(-) diff --git a/src/SIM/SIMoutput.C b/src/SIM/SIMoutput.C index 4f27b7142..a0ad47f59 100644 --- a/src/SIM/SIMoutput.C +++ b/src/SIM/SIMoutput.C @@ -360,6 +360,21 @@ bool SIMoutput::parseOutputTag (const tinyxml2::XMLElement* elem) utl::getAttribute(elem,"precision",myPrec); utl::getAttribute(elem,"vtfsize",myPtSize); + const tinyxml2::XMLElement* comps = elem->FirstChildElement("components"); + if (comps && comps->FirstChild()) + { + char* scmp = strdup(comps->FirstChild()->Value()); + for (char* ccmp = strtok(scmp,"|"); ccmp; ccmp = strtok(nullptr,"|")) + wantComps.insert(ccmp); + free(scmp); + if (!wantComps.empty()) + { + IFEM::cout <<"\tSolution components to output:"; + for (const std::string& s : wantComps) IFEM::cout <<" \""<< s <<"\""; + IFEM::cout << std::endl; + } + } + if (std::string fname; utl::getAttribute(elem,"file",fname)) this->setPointResultFile(fname,elem->FirstChildElement("dump_coordinates")); @@ -403,6 +418,13 @@ void SIMoutput::preprocessResultPoints () { for (ResPtPair& rptp : myPoints) this->preprocessResPtGroup(rptp.first,rptp.second); + + if (myProblem && !myPoints.empty()) + { + IFEM::cout <<"\nSecondary solution components (sol2):"; + for (size_t i = 0; i < myProblem->getNoFields(2); i++) + IFEM::cout <<" \""<< myProblem->getField2Name(i) <<"\""; + } } @@ -892,7 +914,7 @@ bool SIMoutput::writeGlvNo (int& nBlock, int& idBlock, bool SIMoutput::writeGlvT (int iStep, int& geoBlk, int& nBlock) const { - if (myVtf && myProblem->hasTractionValues()) + if (myVtf && myProblem && myProblem->hasTractionValues()) { if (msgLevel > 1) IFEM::cout <<"Writing boundary tractions"<< std::endl; @@ -2464,10 +2486,10 @@ bool SIMoutput::saveResults (const Vectors& psol, double time, int step) const IntVec points, elms; Vec3Vec Xp; Matrix sol, sol2; - std::vector compNames; - std::vector* cmpNames = step > 1 ? nullptr : &compNames; + std::vector cmpNames; for (const ASMbase* pch : myModel) - if (!this->evalResults(psol,gridPts,pch,points,elms,Xp,sol,sol2,cmpNames)) + if (!this->evalResults(psol,gridPts,pch,points,elms,Xp,sol,sol2, + step > 1 && wantComps.empty() ? nullptr : &cmpNames)) return false; if (!sol2.empty() && !sol.augmentRows(sol2)) @@ -2481,10 +2503,14 @@ bool SIMoutput::saveResults (const Vectors& psol, double time, int step) const // Write one file for each result component for (size_t r = 1; r <= sol.rows(); r++) { + // Check if only a subset of the result components are wanted + if (!wantComps.empty() && r <= cmpNames.size() && + wantComps.find(cmpNames[r-1]) == wantComps.end()) continue; + std::string fName = fileName + "_" + std::to_string(r) + fileExt; std::ofstream fs(fName,step > 1 ? std::ios::app : std::ios::out); - if (step == 1 && r <= compNames.size()) - fs <<"# "<< compNames[r-1] << '\n'; + if (step == 1 && r <= cmpNames.size()) + fs <<"# "<< cmpNames[r-1] << '\n'; fs << time; for (size_t c = 1; c <= sol.cols(); c++) fs <<" "<< sol(r,c); diff --git a/src/SIM/SIMoutput.h b/src/SIM/SIMoutput.h index 684ae55c6..37a293c45 100644 --- a/src/SIM/SIMoutput.h +++ b/src/SIM/SIMoutput.h @@ -428,6 +428,8 @@ class SIMoutput : public SIMinput std::map myAddScalars; //!< Scalar functions to output + std::set wantComps; //!< Component names for grid output + int myPrec; //!< Output precision for result sampling double myPtSize; //!< Size of result point visualization in VTF-file int myGeomID; //!< Geometry block ID for the first patch in the VTF-file