From 21a78ada4296c1453dac5eab849ead379bd2e334 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Wed, 5 Nov 2025 09:03:54 +0100 Subject: [PATCH 01/11] Added: ASMLagBase::getBoundingBox(). Changed: Use RealArray instead of Vector when no algebraic operations. --- src/ASM/ASMLagBase.C | 25 +++++++++++++++++++++---- src/ASM/ASMLagBase.h | 6 ++++-- 2 files changed, 25 insertions(+), 6 deletions(-) diff --git a/src/ASM/ASMLagBase.C b/src/ASM/ASMLagBase.C index f38ef3c46..0c8264352 100644 --- a/src/ASM/ASMLagBase.C +++ b/src/ASM/ASMLagBase.C @@ -12,6 +12,7 @@ //============================================================================== #include "ASMLagBase.h" +#include "Vec3Oper.h" ASMLagBase::ASMLagBase (const ASMLagBase& patch, bool cpCoord) : coord(myCoord) @@ -20,7 +21,8 @@ ASMLagBase::ASMLagBase (const ASMLagBase& patch, bool cpCoord) : coord(myCoord) } -bool ASMLagBase::nodalField (Matrix& field, const Vector& sol, size_t nno) const +bool ASMLagBase::nodalField (Matrix& field, + const RealArray& sol, size_t nno) const { size_t nPnts = coord.size(); size_t nComp = sol.size() / nno; @@ -28,7 +30,7 @@ bool ASMLagBase::nodalField (Matrix& field, const Vector& sol, size_t nno) const return false; field.resize(nComp,nPnts); - const double* u = sol.ptr(); + const double* u = sol.data(); for (size_t iPt = 1; iPt <= nPnts; iPt++, u += nComp) field.fillColumn(iPt,u); @@ -48,7 +50,22 @@ Vec3 ASMLagBase::getGeometricCenter (const std::vector& MNPC) const } -bool ASMLagBase::updateCoords (const Vector& displ, unsigned char nsd) +double ASMLagBase::getBoundingBox (const std::vector& MNPC, + Vec3Pair& bbox) const +{ + bbox.first = bbox.second = coord[MNPC.front()]; + for (size_t n = 1; n < MNPC.size(); n++) + for (int i = 0; i < 3; i++) + if (double x = coord[MNPC[n]][i]; x < bbox.first[i]) + bbox.first[i] = x; + else if (x > bbox.second[i]) + bbox.second[i] = x; + + return (bbox.second - bbox.first).length(); +} + + +bool ASMLagBase::updateCoords (const RealArray& displ, unsigned char nsd) { if (displ.size() != nsd*coord.size()) { @@ -58,7 +75,7 @@ bool ASMLagBase::updateCoords (const Vector& displ, unsigned char nsd) return false; } - const double* dpt = displ.ptr(); + const double* dpt = displ.data(); for (Vec3& XYZ : myCoord) { XYZ += RealArray(dpt,dpt+nsd); diff --git a/src/ASM/ASMLagBase.h b/src/ASM/ASMLagBase.h index dc991cb09..db0424113 100644 --- a/src/ASM/ASMLagBase.h +++ b/src/ASM/ASMLagBase.h @@ -34,15 +34,17 @@ class ASMLagBase virtual ~ASMLagBase() {} //! \brief Direct nodal evaluation of a solution field. - bool nodalField(Matrix& field, const Vector& sol, size_t nno) const; + bool nodalField(Matrix& field, const RealArray& sol, size_t nno) const; //! \brief Returns the geometric center of an element. Vec3 getGeometricCenter(const std::vector& MNPC) const; + //! \brief Returns the bounding box and diameter of an element. + double getBoundingBox(const std::vector& MNPC, Vec3Pair& bbox) const; //! \brief Updates the nodal coordinates for this patch. //! \param[in] displ Incremental displacements to update the coordinates with //! \param[in] nsd Number of space dimensions - bool updateCoords(const Vector& displ, unsigned char nsd); + bool updateCoords(const RealArray& displ, unsigned char nsd); public: //! \brief Updates patch origin by adding a constant to all nodes. From 57817471b81f719de12f450057c1b60c445fac79 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Wed, 5 Nov 2025 09:18:05 +0100 Subject: [PATCH 02/11] 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 b6550835007f68d86b032657a8c50150a097374c Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Wed, 5 Nov 2025 10:46:29 +0100 Subject: [PATCH 03/11] Changed: Let the ASMu2DLagTest constructor setup and read the xml input to reduce repetive code. The read() call will not fail anyway (and if it does, it will be caugth by the subsequent checks). Also merge methods genThreadGroups() and getThreadGroups(). Use loops instead of lambdas in the tests for simplicity. Fixed: Swap node 3 and 4 when reading 4-noded element topology. --- src/ASM/Test/TestASMu2DLag.C | 260 ++++++++++++++++------------------- 1 file changed, 115 insertions(+), 145 deletions(-) diff --git a/src/ASM/Test/TestASMu2DLag.C b/src/ASM/Test/TestASMu2DLag.C index 672b2d8a8..c045e8ae9 100644 --- a/src/ASM/Test/TestASMu2DLag.C +++ b/src/ASM/Test/TestASMu2DLag.C @@ -10,184 +10,162 @@ //! //============================================================================== -#include "ASMu2DLag.h" - #include "Catch2Support.h" +#include "ASMu2DLag.h" #include "MPC.h" #include -namespace -{ -class ASMu2DLagTest : public ASMu2DLag +namespace { -public: - ASMu2DLagTest() : ASMu2DLag(2,2,'x') { ASMbase::resetNumbering(); } - - void genThreadGroups(bool separateGroup1noded = false) + std::istream& generateXMLModel(std::stringstream& str, + float Lx, float Ly, int nx, int ny, + int oneNodeElms = 0) { - this->generateThreadGroupsMultiColored(true, separateGroup1noded); - } - const ThreadGroups& getThreadGroups() const { return threadGroups; } -}; - - -void generateXMLModel(std::stringstream& str, - float Lx, float Ly, int nx, int ny, int oneNodeElms = 0) -{ - str << "\n\n"; - const double dx = Lx / nx; - const double dy = Ly / ny; - for (int j = 0; j < ny+1; ++j) - for (int i = 0; i < nx+1; ++i) - str << i*dx << " " << j*dy << " 0.0\n"; - for (int i = 0; i < oneNodeElms; ++i) - str << i+2 << " 0 0\n"; - str << "\n\n"; - for (int j = 0; j < ny; ++j) - for (int i = 0; i < nx; ++i) - str << i + j*(nx+1) << " " << i + 1 + j*(nx+1) << " " - << i + 1 + (j+1)*(nx+1) << " " << i + (j+1)*(nx+1) << '\n'; - str << "\n"; - if (oneNodeElms > 0) { - str << "\n"; + str << "\n\n"; + const double dx = Lx / nx; + const double dy = Ly / ny; + for (int j = 0; j <= ny; ++j) + for (int i = 0; i <= nx; ++i) + str << i*dx << " " << j*dy << " 0.0\n"; for (int i = 0; i < oneNodeElms; ++i) - str << i+(nx+1)*(ny+1) << '\n'; + str << i+2 << " 0 0\n"; + str << "\n\n"; + for (int j = 0; j < ny; ++j) + for (int i = 0; i < nx; ++i) + str << i + j *(nx+1) << " " << i+1 + j *(nx+1) << " " + << i + (j+1)*(nx+1) << " " << i+1 + (j+1)*(nx+1) << "\n"; str << "\n"; + if (oneNodeElms > 0) { + str << "\n"; + for (int i = 0; i < oneNodeElms; ++i) + str << (nx+1)*(ny+1) + i << "\n"; + str << "\n"; + } + str << "\n"; + return str; } - str << "\n"; -} + class ASMu2DLagTest : public ASMu2DLag + { + public: + ASMu2DLagTest(float Lx, float Ly, int nx, int ny, int oneNodeElms = 0) + : ASMu2DLag(2,2,'x') + { + ASMbase::resetNumbering(); + std::stringstream myStream; + this->read(generateXMLModel(myStream,Lx,Ly,nx,ny,oneNodeElms)); + } + + const IntMat& genThreadGroups(bool separateGroup1noded = false) + { + this->generateThreadGroupsMultiColored(false,separateGroup1noded); + return threadGroups[0]; + } + }; } TEST_CASE("TestASMu2DLag.GenerateThreadGroups3x3") { - std::stringstream str; - generateXMLModel(str, 1.0, 1.0, 3, 3); - - ASMu2DLagTest pch; - REQUIRE(pch.read(str)); + ASMu2DLagTest pch(1.0, 1.0, 3, 3); REQUIRE(pch.generateFEMTopology()); - pch.genThreadGroups(); - const ThreadGroups& groups = pch.getThreadGroups(); - REQUIRE(groups[0].size() == 4); - REQUIRE(groups[0][0].size() == 4); - REQUIRE(groups[0][1].size() == 2); - REQUIRE(groups[0][2].size() == 2); - REQUIRE(groups[0][3].size() == 1); - REQUIRE(groups[0][0] == std::vector{0, 2, 6, 8}); - REQUIRE(groups[0][1] == std::vector{1, 7}); - REQUIRE(groups[0][2] == std::vector{3, 5,}); - REQUIRE(groups[0][3] == std::vector{4}); + const IntMat& group = pch.genThreadGroups(); + + REQUIRE(group.size() == 4); + REQUIRE(group[0].size() == 4); + REQUIRE(group[1].size() == 2); + REQUIRE(group[2].size() == 2); + REQUIRE(group[3].size() == 1); + REQUIRE(group[0] == std::vector{0, 2, 6, 8}); + REQUIRE(group[1] == std::vector{1, 7}); + REQUIRE(group[2] == std::vector{3, 5,}); + REQUIRE(group[3] == std::vector{4}); } TEST_CASE("TestASMu2DLag.GenerateThreadGroups4x4") { - std::stringstream str; - generateXMLModel(str, 1.0, 1.0, 4, 4); - - ASMu2DLagTest pch; - REQUIRE(pch.read(str)); + ASMu2DLagTest pch(1.0, 1.0, 4, 4); REQUIRE(pch.generateFEMTopology()); - pch.genThreadGroups(); - const ThreadGroups& groups = pch.getThreadGroups(); - REQUIRE(groups[0].size() == 4); - REQUIRE(groups[0][0].size() == 4); - REQUIRE(groups[0][1].size() == 4); - REQUIRE(groups[0][2].size() == 4); - REQUIRE(groups[0][3].size() == 4); - REQUIRE(groups[0][0] == std::vector{0, 2, 8, 10}); - REQUIRE(groups[0][1] == std::vector{1, 3, 9, 11}); - REQUIRE(groups[0][2] == std::vector{4, 6, 12, 14}); - REQUIRE(groups[0][3] == std::vector{5, 7, 13, 15}); + const IntMat& group = pch.genThreadGroups(); + + REQUIRE(group.size() == 4); + REQUIRE(group[0].size() == 4); + REQUIRE(group[1].size() == 4); + REQUIRE(group[2].size() == 4); + REQUIRE(group[3].size() == 4); + REQUIRE(group[0] == std::vector{0, 2, 8, 10}); + REQUIRE(group[1] == std::vector{1, 3, 9, 11}); + REQUIRE(group[2] == std::vector{4, 6, 12, 14}); + REQUIRE(group[3] == std::vector{5, 7, 13, 15}); } TEST_CASE("TestASMu2DLag.GenerateThreadGroups3x3TwoPC") { - std::stringstream str; - generateXMLModel(str, 1.0, 1.0, 3, 3); - - ASMu2DLagTest pch; - REQUIRE(pch.read(str)); - pch.add2PC(1, 1, 16); - REQUIRE(pch.generateFEMTopology()); + ASMu2DLagTest pch1(1.0, 1.0, 3, 3); + pch1.add2PC(1, 1, 16); + REQUIRE(pch1.generateFEMTopology()); - str.clear(); - str.seekg(0, std::ios_base::beg); - ASMu2DLagTest pch2; - REQUIRE(pch2.read(str)); + ASMu2DLagTest pch2(1.0, 1.0, 3, 3); pch2.add2PC(16, 1, 1); REQUIRE(pch2.generateFEMTopology()); - auto checks = [](ASMu2DLagTest& p) + for (ASMu2DLagTest* pch : { &pch1, &pch2 }) { - p.genThreadGroups(); - const ThreadGroups& groups = p.getThreadGroups(); - REQUIRE(groups[0].size() == 5); - REQUIRE(groups[0][0].size() == 3); - REQUIRE(groups[0][1].size() == 2); - REQUIRE(groups[0][2].size() == 2); - REQUIRE(groups[0][3].size() == 1); - REQUIRE(groups[0][4].size() == 1); - REQUIRE(groups[0][0] == std::vector{0, 2, 6}); - REQUIRE(groups[0][1] == std::vector{1, 7}); - REQUIRE(groups[0][2] == std::vector{3, 5,}); - REQUIRE(groups[0][3] == std::vector{4}); - REQUIRE(groups[0][4] == std::vector{8}); - }; - - checks(pch); - checks(pch2); + const IntMat& group = pch->genThreadGroups(); + + REQUIRE(group.size() == 5); + REQUIRE(group[0].size() == 3); + REQUIRE(group[1].size() == 2); + REQUIRE(group[2].size() == 2); + REQUIRE(group[3].size() == 1); + REQUIRE(group[4].size() == 1); + REQUIRE(group[0] == std::vector{0, 2, 6}); + REQUIRE(group[1] == std::vector{1, 7}); + REQUIRE(group[2] == std::vector{3, 5,}); + REQUIRE(group[3] == std::vector{4}); + REQUIRE(group[4] == std::vector{8}); + } } TEST_CASE("TestASMu2DLag.GenerateThreadGroups3x3MPC") { - std::stringstream str; - generateXMLModel(str, 1.0, 1.0, 3, 3); - - ASMu2DLagTest pch; - REQUIRE(pch.read(str)); + ASMu2DLagTest pch(1.0, 1.0, 3, 3); MPC* mpc = new MPC(1, 1); mpc->addMaster(16, 1); mpc->addMaster(13, 1); pch.addMPC(mpc); REQUIRE(pch.generateFEMTopology()); - pch.genThreadGroups(); - - const ThreadGroups& groups = pch.getThreadGroups(); - REQUIRE(groups[0].size() == 4); - REQUIRE(groups[0][0].size() == 3); - REQUIRE(groups[0][1].size() == 3); - REQUIRE(groups[0][2].size() == 2); - REQUIRE(groups[0][3].size() == 1); - REQUIRE(groups[0][0] == std::vector{0, 2, 7}); - REQUIRE(groups[0][1] == std::vector{1, 6, 8}); - REQUIRE(groups[0][2] == std::vector{3, 5,}); - REQUIRE(groups[0][3] == std::vector{4}); + const IntMat& group = pch.genThreadGroups(); + + REQUIRE(group.size() == 4); + REQUIRE(group[0].size() == 3); + REQUIRE(group[1].size() == 3); + REQUIRE(group[2].size() == 2); + REQUIRE(group[3].size() == 1); + REQUIRE(group[0] == std::vector{0, 2, 7}); + REQUIRE(group[1] == std::vector{1, 6, 8}); + REQUIRE(group[2] == std::vector{3, 5,}); + REQUIRE(group[3] == std::vector{4}); } TEST_CASE("TestASMu2DLag.GenerateThreadGroups3x3OneNode") { - std::stringstream str; - generateXMLModel(str, 1.0, 1.0, 3, 3, 4); - - ASMu2DLagTest pch; - REQUIRE(pch.read(str)); + ASMu2DLagTest pch(1.0, 1.0, 3, 3, 4); REQUIRE(pch.generateFEMTopology()); - auto checks = [](ASMu2DLagTest& p, bool with1) + for (bool with1 : { false, true }) { - p.genThreadGroups(with1); - const ThreadGroups& groups = p.getThreadGroups(); + const IntMat& group = pch.genThreadGroups(with1); + const auto ref = with1 ? std::vector{ std::vector{9, 10, 11, 12}, @@ -203,32 +181,26 @@ TEST_CASE("TestASMu2DLag.GenerateThreadGroups3x3OneNode") std::vector{3, 5}, std::vector{4} }; - REQUIRE(groups[0].size() == ref.size()); + + REQUIRE(group.size() == ref.size()); for (size_t i = 0; i < ref.size(); ++i) { - REQUIRE(groups[0][i].size() == ref[i].size()); - REQUIRE(groups[0][i] == ref[i]); + REQUIRE(group[i].size() == ref[i].size()); + REQUIRE(group[i] == ref[i]); } - }; - - checks(pch, false); - checks(pch, true); + } } TEST_CASE("TestASMu2DLag.GenerateThreadGroups3x3OneNodeSPC") { - std::stringstream str; - generateXMLModel(str, 1.0, 1.0, 3, 3, 4); - - ASMu2DLagTest pch; - REQUIRE(pch.read(str)); + ASMu2DLagTest pch(1.0, 1.0, 3, 3, 4); pch.add2PC(17, 1, 1); REQUIRE(pch.generateFEMTopology()); - auto checks = [](ASMu2DLagTest& p, bool with1) + for (bool with1 : { false, true }) { - p.genThreadGroups(with1); - const ThreadGroups& groups = p.getThreadGroups(); + const IntMat& group = pch.genThreadGroups(with1); + const auto ref = with1 ? std::vector{ std::vector{10, 11, 12}, @@ -244,13 +216,11 @@ TEST_CASE("TestASMu2DLag.GenerateThreadGroups3x3OneNodeSPC") std::vector{3, 5}, std::vector{4} }; - REQUIRE(groups[0].size() == ref.size()); + + REQUIRE(group.size() == ref.size()); for (size_t i = 0; i < ref.size(); ++i) { - REQUIRE(groups[0][i].size() == ref[i].size()); - REQUIRE(groups[0][i] == ref[i]); + REQUIRE(group[i].size() == ref[i].size()); + REQUIRE(group[i] == ref[i]); } - }; - - checks(pch, false); - checks(pch, true); + } } From 7c46a23d3cd6e603c31c99e051dfeda67c1caba2 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Thu, 6 Nov 2025 08:57:26 +0100 Subject: [PATCH 04/11] Added: Fortran code for parametric inversion of 3- and 4-noded shells --- CMakeLists.txt | 13 +- src/ASM/Ftn/dcros3.f | 12 ++ src/ASM/Ftn/dmuax3.f | 12 ++ src/ASM/Ftn/dnorm3.f | 33 +++++ src/ASM/Ftn/dprint.f | 49 +++++++ src/ASM/{ => Ftn}/dslbln.f | 0 src/ASM/{ => Ftn}/dsolv3.f | 0 src/ASM/Ftn/errmsg.f | 42 ++++++ src/ASM/Ftn/q4_inv.f | 274 +++++++++++++++++++++++++++++++++++++ src/ASM/Ftn/shape2.f | 34 +++++ src/ASM/Ftn/shltri.f | 190 +++++++++++++++++++++++++ src/ASM/Ftn/t3_inv.f | 216 +++++++++++++++++++++++++++++ 12 files changed, 873 insertions(+), 2 deletions(-) create mode 100644 src/ASM/Ftn/dcros3.f create mode 100644 src/ASM/Ftn/dmuax3.f create mode 100644 src/ASM/Ftn/dnorm3.f create mode 100644 src/ASM/Ftn/dprint.f rename src/ASM/{ => Ftn}/dslbln.f (100%) rename src/ASM/{ => Ftn}/dsolv3.f (100%) create mode 100644 src/ASM/Ftn/errmsg.f create mode 100644 src/ASM/Ftn/q4_inv.f create mode 100644 src/ASM/Ftn/shape2.f create mode 100644 src/ASM/Ftn/shltri.f create mode 100644 src/ASM/Ftn/t3_inv.f diff --git a/CMakeLists.txt b/CMakeLists.txt index abf3ee350..182ad0939 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -100,8 +100,17 @@ target_sources( src/ASM/BDFMats.C src/ASM/BlockElmMats.C src/ASM/DomainDecomposition.C - src/ASM/dslbln.f - src/ASM/dsolv3.f + src/ASM/Ftn/dcros3.f + src/ASM/Ftn/dmuax3.f + src/ASM/Ftn/dnorm3.f + src/ASM/Ftn/dprint.f + src/ASM/Ftn/dslbln.f + src/ASM/Ftn/dsolv3.f + src/ASM/Ftn/errmsg.f + src/ASM/Ftn/q4_inv.f + src/ASM/Ftn/shape2.f + src/ASM/Ftn/shltri.f + src/ASM/Ftn/t3_inv.f src/ASM/DualField.C src/ASM/ElmMats.C src/ASM/Field.C diff --git a/src/ASM/Ftn/dcros3.f b/src/ASM/Ftn/dcros3.f new file mode 100644 index 000000000..2cc69e287 --- /dev/null +++ b/src/ASM/Ftn/dcros3.f @@ -0,0 +1,12 @@ + SUBROUTINE DCROS3 (X,Y,Z) +C +C Calculate {Z} = {X} x {Y} +C + DOUBLE PRECISION X(3),Y(3),Z(3) +C + Z(1) = X(2)*Y(3) - Y(2)*X(3) + Z(2) = X(3)*Y(1) - Y(3)*X(1) + Z(3) = X(1)*Y(2) - Y(1)*X(2) +C + RETURN + END diff --git a/src/ASM/Ftn/dmuax3.f b/src/ASM/Ftn/dmuax3.f new file mode 100644 index 000000000..582f6f1fc --- /dev/null +++ b/src/ASM/Ftn/dmuax3.f @@ -0,0 +1,12 @@ + SUBROUTINE DMUAX3 (A,X,Y) +C +C Calculate {Y} = [A] * {X} +C + DOUBLE PRECISION A(3,3),X(3),Y(3) +C + Y(1) = A(1,1)*X(1) + A(1,2)*X(2) + A(1,3)*X(3) + Y(2) = A(2,1)*X(1) + A(2,2)*X(2) + A(2,3)*X(3) + Y(3) = A(3,1)*X(1) + A(3,2)*X(2) + A(3,3)*X(3) +C + RETURN + END diff --git a/src/ASM/Ftn/dnorm3.f b/src/ASM/Ftn/dnorm3.f new file mode 100644 index 000000000..441969c69 --- /dev/null +++ b/src/ASM/Ftn/dnorm3.f @@ -0,0 +1,33 @@ + DOUBLE PRECISION FUNCTION DNORM3 (V) + DOUBLE PRECISION V(3),DV +C +C DNORM3 : Calculate |v| and 1/|v| * {v} +C + DOUBLE PRECISION ZERO , EPS , ONE + PARAMETER ( ZERO = 0.0D0, EPS = 1.0D-8, ONE = 1.0D0 ) +C + DV = SQRT(V(1)*V(1) + V(2)*V(2) + V(3)*V(3)) + IF (DV .LT. EPS) THEN + DNORM3 = ZERO + ELSE + DNORM3 = DV + V(1) = V(1) / DV + V(2) = V(2) / DV + V(3) = V(3) / DV + IF (ABS(V(1)) .GT. ONE-EPS) THEN + V(1) = SIGN(ONE,V(1)) + V(2) = ZERO + V(3) = ZERO + ELSE IF (ABS(V(2)) .GT. ONE-EPS) THEN + V(1) = ZERO + V(2) = SIGN(ONE,V(2)) + V(3) = ZERO + ELSE IF (ABS(V(3)) .GT. ONE-EPS) THEN + V(1) = ZERO + V(2) = ZERO + V(3) = SIGN(ONE,V(3)) + ENDIF + ENDIF +C + RETURN + END diff --git a/src/ASM/Ftn/dprint.f b/src/ASM/Ftn/dprint.f new file mode 100644 index 000000000..383baae62 --- /dev/null +++ b/src/ASM/Ftn/dprint.f @@ -0,0 +1,49 @@ + SUBROUTINE DPRINT (A,M,N,CHEAD,IWR) + INTEGER M,N, IWR + DOUBLE PRECISION A(M,*) + CHARACTER CHEAD*(*) +C +C PURPOSE: +C DPRINT : Print the double precision array A(M,N). +C +C METHOD: +C The array is printed as a two-dimensional matrix. +C Rows and columns are numbered and each element is printed with +C the format 1PE14.6, i.e. 7 significant digits. +C A heading (e.g. the array name) is also printed. +C +C INPUT ARGUMENTS: +C A - Array to be printed +C M - Number of rows +C N - Number of columns +C CHEAD - Heading +C IWR - Print unit number +C + INTEGER NCOL,J,I,K,L + PARAMETER ( NCOL = 9 ) +C + IF (N .LE. 0 .OR. M .LE. 0) GOTO 8000 +C + WRITE(IWR,9000) CHEAD +C + DO 2000 J = 1,N,NCOL + K = MIN(N,J+NCOL-1) +C + WRITE(IWR,9100) (L,L=J,K) + WRITE(IWR,9200) ('--------------',L=J,K) +C + DO 1000 I = 1,M + WRITE(IWR,9300) I,(A(I,L),L=J,K) + 1000 CONTINUE +C + 2000 CONTINUE +C + 8000 CONTINUE + RETURN +C + 9000 FORMAT( / 5X, A, ' = ' ) + 9100 FORMAT( / 4X, 9(I9,5X) ) + 9200 FORMAT( 4X, 9A14 ) + 9300 FORMAT( I3, '|', 1P,9E14.6 ) +C + END diff --git a/src/ASM/dslbln.f b/src/ASM/Ftn/dslbln.f similarity index 100% rename from src/ASM/dslbln.f rename to src/ASM/Ftn/dslbln.f diff --git a/src/ASM/dsolv3.f b/src/ASM/Ftn/dsolv3.f similarity index 100% rename from src/ASM/dsolv3.f rename to src/ASM/Ftn/dsolv3.f diff --git a/src/ASM/Ftn/errmsg.f b/src/ASM/Ftn/errmsg.f new file mode 100644 index 000000000..9851b0dea --- /dev/null +++ b/src/ASM/Ftn/errmsg.f @@ -0,0 +1,42 @@ + SUBROUTINE ERRMSG (PRGNAM,IWR,IERR) + CHARACTER PRGNAM*(*) + INTEGER IWR,IERR +C +C ERRMSG : Print an error/warning message. +C +C PRGNAM - Name of subroutine where the error occurred +C IERR - Error flag +C > 1 : Warning +C = -1 : Traceback +C < -1 : Error +C + IF (IERR .LT. -1) THEN + WRITE(IWR,9901) PRGNAM + ELSE IF (IERR .GT. 1) THEN + WRITE(IWR,9902) PRGNAM + ELSE IF (IERR .EQ. -1) THEN + WRITE(IWR,9903) PRGNAM + ENDIF +C + IF (ABS(IERR) .LE. 1) GOTO 8000 +C + IF (IERR .EQ. -2) THEN + WRITE(IWR,9904) 'INVALID INPUT ARGUMENTS' + ELSE IF (IERR .EQ. -98) THEN + WRITE(IWR,9904) 'ARRAY IS NOT ALLOCATED' + ELSE IF (IERR .EQ. -99) THEN + WRITE(IWR,9904) 'UNABLE TO ALLOCATE MEMORY' + ELSE + WRITE(IWR,9905) IERR + ENDIF +C + 8000 CONTINUE + RETURN +C + 9901 FORMAT( / ' *** ERROR RETURN FROM SUBROUTINE ', A ) + 9902 FORMAT( / ' ** WARNING FROM SUBROUTINE ', A ) + 9903 FORMAT( ' * TRACEBACK: ', A ) + 9904 FORMAT( 5X, A ) + 9905 FORMAT( 5X, 'ERROR FLAG = ', I6 ) +C + END diff --git a/src/ASM/Ftn/q4_inv.f b/src/ASM/Ftn/q4_inv.f new file mode 100644 index 000000000..77dcb674c --- /dev/null +++ b/src/ASM/Ftn/q4_inv.f @@ -0,0 +1,274 @@ + SUBROUTINE Q4_INV (NDIM,EPSO,EPSZ,X,TENC,XI,IPSW,IWR,IERR) +C +C $Id$ +C***************************** USERS' PART ***************************** +C +C A F E M ROUTINE : Q4_INV +C +C PURPOSE: +C Q4_INV : Find natural coordinates for a point within a Q4-element. +C +C METHOD: +C The natural coordinates {XI} corresponding to the +C point {X} are calculated. The element nodes are assumed to have +C natural coordinates (-1,-1), (1,-1), (1,1), (-1,1), respectively. +C +C INPUT ARGUMENTS: +C NDIM - Number of spatial dimensions +C EPSO - Tolerence for checking if point is outside the element +C EPSZ - Tolerance for equality check with zero +C X - Cartesian coordinates of a point within the element +C TENC - Table of element nodal coordinates +C IPSW - Print switch +C IWR - Print unit number +C +C INPUT/OUTPUT ARGUMENTS: +C None +C +C OUTPUT ARGUMENTS: +C XI - Natural coordinates of the point {X} +C IERR - Error flag +C +C RESTRICTIONS: +C None +C +C************************** MAINTENANCE PART *************************** +C +C INTERNAL VARIABLES: +C None +C +C PRINT SWITCH: +C .LT. 0 - No print +C .EQ. 0 - Print error messages only +C .GE. 7 - Print ENTER and LEAVE +C .GE. 8 - Print in addition INPUT parameters +C .GE. 9 - Print in addition OUTPUT parameters +C +C SUBROUTINES CALLED: +C DDOT (double precision function) +C DMUAX3 +C DSLBLN +C SHAPE2 +C SHLTRI +C +C CODING: +C 19-SEP-96, Knut Morten Okstad, NTNU +C +C REVISION: +C +C********************** END SUBROUTINE DESCRIPTION ********************* +C + IMPLICIT NONE +C +C * Global variables +C + INTEGER NDIM,IPSW,IWR,IERR + DOUBLE PRECISION EPSO,EPSZ,X(NDIM),TENC(NDIM,4),XI(2) +C +C * Internal variables + LOGICAL LXY + INTEGER i,IERRL,j,NSOL + DOUBLE PRECISION A(4,2),D1,D2,TOL,TGL(9),TLX(2,5), + + X1(4),X2(4),XL(3) +C + DOUBLE PRECISION ZERO , C , ONE + PARAMETER ( ZERO = 0.0D0, C = 0.25D0, ONE = 1.0D0 ) +C + DOUBLE PRECISION DDOT +C +C**************************** END VARIABLES **************************** +C +C Entry section +C + IF (IPSW .GE. 7) WRITE(IWR,9900) + IF (IPSW .GE. 8) THEN + WRITE(IWR,9920) EPSO,EPSZ,X + CALL DPRINT (TENC,NDIM,4,'TENC ',IWR) + ENDIF +C +C +C Program logic section +C +C * Check input arguments + IERR = -1 + IF (EPSZ .LE. ZERO .OR. EPSO .LT. EPSZ) GOTO 7200 +C + IF (NDIM .EQ. 2) THEN + LXY = .TRUE. + ELSE IF (NDIM .EQ. 3) THEN + TOL = MAX(TENC(3,1),TENC(3,2),TENC(3,3),TENC(3,4)) + + - MIN(TENC(3,1),TENC(3,2),TENC(3,3),TENC(3,4)) + LXY = TOL .LT. EPSZ + ELSE + GOTO 7200 + ENDIF +C + IF (LXY) THEN +C +C * The element is parallel +C to the global XY-plane + DO 1000 j = 1,2 + A(1,j) = ( TENC(j,1) - TENC(j,2) + TENC(j,3) - TENC(j,4))*C + A(2,j) = (-TENC(j,1) + TENC(j,2) + TENC(j,3) - TENC(j,4))*C + A(3,j) = (-TENC(j,1) - TENC(j,2) + TENC(j,3) + TENC(j,4))*C + A(4,j) = ( TENC(j,1) + TENC(j,2) + TENC(j,3) + TENC(j,4))*C + 1000 CONTINUE +C + A(4,1) = X(1) - A(4,1) + A(4,2) = X(2) - A(4,2) +C + ELSE +C +C * The element is not parallel +C to the global XY-plane, need +C to calculate local coordinates + DO 2000 j = 1,NDIM + X1(j) = TENC(j,3) - TENC(j,1) + X2(j) = TENC(j,4) - TENC(j,2) + 2000 CONTINUE + CALL SHLTRI (2,X1,X2,TGL,IPSW,IWR,IERRL) + IF (IERRL .LT. 0) GOTO 7100 +C + DO 2100 j = 1,4 + CALL DMUAX3 (TGL,TENC(1,j),TLX(1,j)) + 2100 CONTINUE + CALL DMUAX3 (TGL,X,XL) +C +C * Initiate the aux. variables + DO 2200 j = 1,2 + A(1,j) = ( TLX(j,1) - TLX(j,2) + TLX(j,3) - TLX(j,4)) * C + A(2,j) = (-TLX(j,1) + TLX(j,2) + TLX(j,3) - TLX(j,4)) * C + A(3,j) = (-TLX(j,1) - TLX(j,2) + TLX(j,3) + TLX(j,4)) * C + A(4,j) = ( TLX(j,1) + TLX(j,2) + TLX(j,3) + TLX(j,4)) * C + 2200 CONTINUE +C + A(4,1) = XL(1) - A(4,1) + A(4,2) = XL(2) - A(4,2) +C + ENDIF +C +C----------------------------------------------------------------------- +C +C We have to solve the following set of equations (j=1,2): +C +C A1j*XI*ETA + A2j*XI + A3j*ETA = A4j +C +C The way that we may solve this nonlinear (in XI and ETA) +C set of equations depends on the coefficients Aij. +C The solution is unique for proper input. See DSLBLN for details. +C +C----------------------------------------------------------------------- +C + CALL DSLBLN (IPSW,IWR,EPSZ,A(1,1),A(1,2),NSOL,X1,X2) + IF (NSOL .LE. 0) GOTO 7100 +C +C * Check that the solutions +C are "inside" the element + TOL = ONE + EPSO + DO 3000 j = 1,NSOL + IF (ABS(X1(j)).LT.TOL .AND. ABS(X2(j)).LT.TOL) THEN + IERR = IERR + 1 + IF (IERR .GT. 0) THEN +C +C * We have multiple solutions. +C Check if they are almost equal. +C + IF (ABS(XI(1)-X1(j)) .LT. EPSZ*(TOL+TOL) .AND. + + ABS(XI(2)-X2(j)) .LT. EPSZ*(TOL+TOL)) THEN + XI(1) = 0.5D0*(XI(1)+X1(j)) + XI(2) = 0.5D0*(XI(2)+X2(j)) + IERR = IERR - 1 + GOTO 3000 + ENDIF +C +C * Choose the one closest to X +C + CALL SHAPE2 (0,4,XI,TGL,IPSW,IWR,IERRL) + IF (IERRL .LT. 0) GOTO 7100 +C + D1 = ZERO + DO 2500 i = 1,NDIM + D1 = D1 + (X(i)-DDOT(4,TGL,1,TENC(i,1),NDIM))**2.0D0 + 2500 CONTINUE +C + XL(1) = X1(j) + XL(2) = X2(j) + CALL SHAPE2 (0,4,XL,TGL,IPSW,IWR,IERRL) + IF (IERRL .LT. 0) GOTO 7100 +C + D2 = ZERO + DO 2600 i = 1,NDIM + D2 = D2 + (X(i)-DDOT(4,TGL,1,TENC(i,1),NDIM))**2.0D0 + 2600 CONTINUE + IF (D1 .LT. D2) GOTO 3000 +C + ENDIF + XI(1) = X1(j) + XI(2) = X2(j) + ENDIF + 3000 CONTINUE + + IF (IERR .GT. 0) GOTO 6200 + IF (IERR .EQ. 0) GOTO 8000 +C +C +C Warning section +C + 6100 CONTINUE + IERR = 1 + IF (IPSW .GE. 0) THEN + WRITE(IWR,9961) (j,X1(j),X2(j),j=1,NSOL) + WRITE(IWR,9920) EPSO,EPSZ,X + CALL DPRINT (TENC,NDIM,4,'TENC ',IWR) + ENDIF + GOTO 8000 + 6200 CONTINUE + IERR = 2 + IF (IPSW .GE. 0) THEN + WRITE(IWR,9962) (j,X1(j),X2(j),j=1,NSOL) + WRITE(IWR,9963) XI,MIN(D1,D2) + WRITE(IWR,9920) EPSO,EPSZ,X + CALL DPRINT (TENC,NDIM,4,'TENC ',IWR) + ENDIF + GOTO 8000 +C +C Error section +C + 7100 CONTINUE + IERR = -1 + GOTO 7900 + 7200 CONTINUE + IERR = -2 +C + 7900 CONTINUE + IF (IPSW .GE. 0) THEN + CALL ERRMSG ('Q4_INV',IWR,IERR) + ENDIF +C +C +C Closing section +C + 8000 CONTINUE + IF (IPSW .GE. 7) WRITE(IWR,9910) + IF (IPSW .GE. 9) WRITE(IWR,9930) IERR,XI +C + RETURN +C + 9900 FORMAT( / ' ENTERING SUBROUTINE Q4_INV' ) + 9910 FORMAT( / ' LEAVING SUBROUTINE Q4_INV' ) + 9920 FORMAT( ' WITH INPUT PARAMETERS:' + + / ' EPS = ', 1P,2E13.5, + + / ' X = ', 3E13.5 ) + 9930 FORMAT( ' WITH OUTPUT PARAMETERS:' + + / ' IERR = ', I6, + + / ' XI,ETA = ', 1P,2E13.5 ) + 9961 FORMAT( / ' ** WARNING FROM SUBROUTINE Q4_INV' + + / ' THE POINT(S) ARE OUTSIDE THE ELEMENT', + + /(' XI(',I1,') = ', 1P,2E13.5) ) + 9962 FORMAT( / ' ** WARNING FROM SUBROUTINE Q4_INV' + + / ' THE SOLUTION IS NOT UNIQUE' + + /(' XI(',I1,') = ', 1P,2E13.5) ) + 9963 FORMAT( ' CHOSE THE SOLUTION CLOSEST TO {X}' + + / ' XI,ETA = ', 1P,2E13.5, ' DELTA = ', E12.5 ) +C + END diff --git a/src/ASM/Ftn/shape2.f b/src/ASM/Ftn/shape2.f new file mode 100644 index 000000000..51907aeed --- /dev/null +++ b/src/ASM/Ftn/shape2.f @@ -0,0 +1,34 @@ + SUBROUTINE SHAPE2 (IOP,NENOD,XI,S,IPSW,IWR,IERR) + IMPLICIT NONE + INTEGER IOP,NENOD,IPSW,IWR,IERR + DOUBLE PRECISION XI(2),S(NENOD) +C + DOUBLE PRECISION UMX1,UMX2,UPX1,UPX2 +C + DOUBLE PRECISION QUART , ONE + PARAMETER ( QUART = 0.25D0, ONE = 1.0D0 ) +C + IF (IOP .EQ. 0 .AND. NENOD .EQ. 3) THEN + IERR = 0 + S(1) = XI(1) + S(2) = XI(2) + S(3) = ONE - XI(1) - XI(2) + ELSE IF (IOP .EQ. 0 .AND. NENOD .EQ. 4) THEN + IERR = 0 + UMX1 = ONE - XI(1) + UMX2 = ONE - XI(2) + UPX1 = ONE + XI(1) + UPX2 = ONE + XI(2) + S(1) = QUART * UMX1*UMX2 + S(2) = QUART * UPX1*UMX2 + S(3) = QUART * UPX1*UPX2 + S(4) = QUART * UMX1*UPX2 + ELSE + IERR = -2 + IF (IPSW .GE. 0) THEN + write(IWR,*) '*** SHAPE2:',IOP,NENOD,'not supported' + ENDIF + ENDIF +C + RETURN + END diff --git a/src/ASM/Ftn/shltri.f b/src/ASM/Ftn/shltri.f new file mode 100644 index 000000000..6b5c8b38c --- /dev/null +++ b/src/ASM/Ftn/shltri.f @@ -0,0 +1,190 @@ + SUBROUTINE SHLTRI (IOP,VX,V3,TGL,IPSW,IWR,IERR) +C +C $Id$ +C***************************** USERS' PART ***************************** +C +C A F E M ROUTINE : SHLTRI +C +C PURPOSE: +C SHLTRI : Calculate a local triad on a shell surface. +C +C METHOD: +C The local Z-axis is defined by the given vector V3. For IOP=0, the +C local X-axis is defined by the projection of the global X-axis +C onto the shell surface. If this projection is zero, the local +C Y-axis is instead defined by the projection of the global Y-axis. +C For IOP=1, the local X-axis is defined by the projection of the +C given vector VX onto the shell surface. For IOP=2, the local +C X-axis is defined by the given vector VX and the local Z-axis is +C defined by the cross product between VX and V3. +C +C INPUT ARGUMENTS: +C IOP - Option telling how to construct the triad (see above) +C VX - Vector defining the local X-axis (for IOP>0 only) +C IPSW - Print switch +C IWR - Print unit number +C +C INPUT/OUTPUT ARGUMENTS: +C V3 - Vector defining the local Z-axis +C +C OUTPUT ARGUMENTS: +C TGL - Transformation matrix +C IERR - Error flag +C +C RESTRICTIONS: +C None +C +C************************** MAINTENANCE PART *************************** +C +C INTERNAL VARIABLES: +C V1, V2 - In-plane base vectors +C +C PRINT SWITCH: +C .LT. 0 - No print +C .EQ. 0 - Print error messages only +C .GE. 6 - Print ENTER and LEAVE +C .GE. 7 - Print in addition INPUT parameters +C .GE. 8 - Print in addition OUTPUT parameters +C +C SUBROUTINES CALLED: +C DCROS3 +C DNORM3 +C +C CODING: +C 3-JUN-96, Knut Morten Okstad, NTNU +C +C REVISION: +C +C********************** END SUBROUTINE DESCRIPTION ********************* +C + IMPLICIT NONE +C +C * Global variables +C + INTEGER IOP,IPSW,IWR,IERR + DOUBLE PRECISION VX(3),V3(3),TGL(3,3) +C +C * Internal variables + INTEGER I,J +C + DOUBLE PRECISION ZERO , EPS + PARAMETER ( ZERO = 0.0D0, EPS = 1.0D-8 ) +C + DOUBLE PRECISION V1(3),V2(3) +C + DOUBLE PRECISION DNORM3 +C +C**************************** END VARIABLES **************************** +C +C Entry section +C + IF (IPSW .GE. 6) WRITE(IWR,9900) + IF (IPSW .GE. 7) THEN + WRITE(IWR,9920) IOP,V3 + IF (IOP .EQ. 1) WRITE(IWR,9921) VX + ENDIF +C +C +C Program logic section +C + i = 1 + IERR = 0 + IF (IOP .NE. 2) THEN + IF (DNORM3(V3) .LE. ZERO) GOTO 7200 + ENDIF +C +C * Set up a local Cartesian basis +C (V1,V2,V3) where V3 coincide +C with the plane normal +C + IF (IOP .EQ. 1) THEN +C +C * Define V1 by the projection of VX +C onto the plane, then V2 = V3 x V1 + CALL DCROS3 (V3,VX,V2) + CALL DCROS3 (V2,V3,V1) +C + ELSE IF (IOP .EQ. 2) THEN +C +C * Define V2 as V1 x V3, +C then V3 = V1 x V2 + CALL DCROS3 (V1,V3,V2) + V3(1) = V2(1) + V3(2) = V2(2) + V3(3) = V2(3) + CALL DCROS3 (V1,V2,V3) + IF (DNORM3(V3) .LE. ZERO) GOTO 7200 +C + ELSE IF (ABS(V3(2)) .GT. EPS .OR. ABS(V3(3)) .GT. EPS) THEN +C +C * Define V1 by projecting the +C global X-axis onto the plane, +C then V2 = V3 x V1 +C + V1(1) = V3(2)*V3(2) + V3(3)*V3(3) + V1(2) = -V3(1)*V3(2) + V1(3) = -V3(1)*V3(3) + CALL DCROS3 (V3,V1,V2) +C + ELSE +C +C * Define V2 by projecting the +C global Y-axis onto the plane, +C then V1 = V2 x V3 + V2(1) = -V3(2)*V3(1) + V2(2) = V3(1)*V3(1) + V3(3)*V3(3) + V2(3) = -V3(2)*V3(3) + CALL DCROS3 (V2,V3,V1) +C + ENDIF +C +C * Normalize the base vectors + i = 2 + IF (DNORM3(V1) .LE. ZERO) GOTO 7200 + IF (DNORM3(V2) .LE. ZERO) GOTO 7200 +C +C * Establish the global-to-local +C transformation matrix, Tgl + DO 1000 I = 1,3 + TGL(1,I) = V1(I) + TGL(2,I) = V2(I) + TGL(3,I) = V3(I) + 1000 CONTINUE +C + GOTO 8000 +C +C Error section +C + 7200 CONTINUE + IERR = -2 + IF (IPSW .GE. 0) THEN + CALL ERRMSG ('SHLTRI',IWR,IERR) + WRITE(IWR,9920) IOP,V3 + IF (IOP .EQ. 1 .OR. IOP .EQ. 2) WRITE(IWR,9921) VX + IF (i .EQ. 2) WRITE(IWR,9922) V1,V2 + ENDIF +C +C +C Closing section +C + 8000 CONTINUE + IF (IPSW .GE. 6) WRITE(IWR,9910) + IF (IPSW .GE. 8) THEN + WRITE(IWR,9930) IERR,((TGL(i,j),j=1,3),i=1,3) + ENDIF +C + RETURN +C + 9900 FORMAT( / ' ENTERING SUBROUTINE SHLTRI' ) + 9910 FORMAT( / ' LEAVING SUBROUTINE SHLTRI' ) + 9920 FORMAT( ' WITH INPUT PARAMETERS:' + + / ' IOP = ', I6, + + / ' V3 = ', 1P,3E13.5 ) + 9921 FORMAT( ' VX = ', 1P,3E13.5 ) + 9922 FORMAT( ' V1 = ', 1P,3E13.5 + + / ' V2 = ', 3E13.5 ) + 9930 FORMAT( ' WITH OUTPUT PARAMETERS:' + + / ' IERR = ', I6, + + / ' TGL = ', 1P,3E13.5, 2(/14X,3E13.5) ) +C + END diff --git a/src/ASM/Ftn/t3_inv.f b/src/ASM/Ftn/t3_inv.f new file mode 100644 index 000000000..24dd6698c --- /dev/null +++ b/src/ASM/Ftn/t3_inv.f @@ -0,0 +1,216 @@ + SUBROUTINE T3_INV (NDIM,EPSO,EPSZ,X,TENC,XI,IPSW,IWR,IERR) +C +C $Id$ +C***************************** USERS' PART ***************************** +C +C A F E M ROUTINE : T3_INV +C +C PURPOSE: +C T3_INV : Find natural coordinates for a point within a triangle. +C +C METHOD: +C The natural (area) coordinates {XI}={A1,A2} corresponding to the +C point {X} are calculated. The element nodes are assumed to have +C the natural coordinates (1,0), (0,1) and (0,0), respectively. +C +C INPUT ARGUMENTS: +C NDIM - Number of spatial dimensions +C EPSO - Tolerence for checking if point is outside the element +C EPSZ - Tolerance for equality check with zero +C X - Cartesian coordinates of a point within the element +C TENC - Table of element nodal coordinates +C IPSW - Print switch +C IWR - Print unit number +C +C INPUT/OUTPUT ARGUMENTS: +C None +C +C OUTPUT ARGUMENTS: +C XI - Natural coordinates of the point {X} +C IERR - Error flag +C +C RESTRICTIONS: +C None +C +C************************** MAINTENANCE PART *************************** +C +C INTERNAL VARIABLES: +C AREA - Element area +C VN - Element normal vector +C V12 - Vector from node 1 to node 2 +C V13 - Vector from node 1 to node 3 +C V1X - Vector from node 1 to the point X +C +C PRINT SWITCH: +C .LT. 0 - No print +C .EQ. 0 - Print error messages only +C .GE. 7 - Print ENTER and LEAVE +C .GE. 8 - Print in addition INPUT parameters +C .GE. 9 - Print in addition OUTPUT parameters +C +C SUBROUTINES CALLED: +C DCROS3 +C DNORM3 +C +C CODING: +C 4-NOV-94, Knut Morten Okstad, NTH +C +C REVISION: +C 2-OCT-96, Knut Morten Okstad, NTNU +C The routine is included in the AFEM package. +C +C********************** END SUBROUTINE DESCRIPTION ********************* +C + IMPLICIT NONE +C +C * Global variables +C + INTEGER NDIM,IPSW,IWR,IERR + DOUBLE PRECISION EPSO,EPSZ,X(NDIM),TENC(NDIM,3),XI(2) +C +C * Internal variables +C + DOUBLE PRECISION AUX(3),A2,A3,AREA,TOL, + + V12(3),V13(3),V1X(3),VN(3) +C + DOUBLE PRECISION ZERO , ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) +C + DOUBLE PRECISION DNORM3 +C +C**************************** END VARIABLES **************************** +C +C Entry section +C + IF (IPSW .GE. 7) WRITE(IWR,9900) + IF (IPSW .GE. 8) THEN + WRITE(IWR,9920) NDIM,EPSO,EPSZ,X + CALL DPRINT (TENC,NDIM,3,'TENC ',IWR) + ENDIF +C +C +C Program logic section +C +C * Check input arguments + IERR = 0 + IF (EPSZ .LE. ZERO .OR. EPSO .LT. EPSZ) GOTO 7200 +C + IF (NDIM .EQ. 2) THEN +C +C * Calculate the element area +C + AREA = ( (TENC(1,1)-TENC(1,3))*(TENC(2,2)-TENC(2,3)) - + + (TENC(1,2)-TENC(1,3))*(TENC(2,1)-TENC(2,3)) ) + IF (AREA .LE. EPSZ) GOTO 7200 +C +C * Calculate the area coordinates +C + XI(1) = ( ( X(1) -TENC(1,3))*(TENC(2,2)-TENC(2,3)) - + + (TENC(1,2)-TENC(1,3))*( X(2) -TENC(2,3)) ) / AREA + XI(2) = ( (TENC(1,1)-TENC(1,3))*( X(2) -TENC(2,3)) - + + ( X(1) -TENC(1,3))*(TENC(2,1)-TENC(2,3)) ) / AREA +C + ELSE IF (NDIM .EQ. 3) THEN +C +C * Calculate the element normal +C vector and the element area + V12(1) = TENC(1,2) - TENC(1,1) + V12(2) = TENC(2,2) - TENC(2,1) + V12(3) = TENC(3,2) - TENC(3,1) + V13(1) = TENC(1,3) - TENC(1,1) + V13(2) = TENC(2,3) - TENC(2,1) + V13(3) = TENC(3,3) - TENC(3,1) + CALL DCROS3 (V12,V13,VN) + AREA = DNORM3(VN) + IF (AREA .LE. EPSZ) GOTO 7200 +C +C * Check the distance from point X +C to the element plane + V1X(1) = X(1) - TENC(1,1) + V1X(2) = X(2) - TENC(2,1) + V1X(3) = X(3) - TENC(3,1) + A3 = VN(1)*V1X(1) + VN(2)*V1X(2) + VN(3)*V1X(3) + IF (ABS(A3) .GT. SQRT(AREA)*EPSO) GOTO 7300 +C +C * Calculate the area coordinates + V1X(1) = V1X(1) - A3*VN(1) + V1X(2) = V1X(2) - A3*VN(2) + V1X(3) = V1X(3) - A3*VN(3) +C + CALL DCROS3 (V12,V1X,AUX) + A3 = SIGN(SQRT(AUX(1)*AUX(1)+AUX(2)*AUX(2)+AUX(3)*AUX(3)), + + AUX(1)*VN(1) +AUX(2)*VN(2) +AUX(3)*VN(3)) + CALL DCROS3 (V1X,V13,AUX) + A2 = SIGN(SQRT(AUX(1)*AUX(1)+AUX(2)*AUX(2)+AUX(3)*AUX(3)), + + AUX(1)*VN(1) +AUX(2)*VN(2) +AUX(3)*VN(3)) +C + XI(1) = (AREA-A2-A3) / AREA + XI(2) = A2 / AREA +C + ELSE + GOTO 7200 + ENDIF +C +C +C * Check that the coordinates are +C "inside" the element + TOL = ONE + EPSO + IF (XI(1) .LT. -EPSO .OR. XI(1) .GT. TOL) GOTO 6100 + IF ( XI(2) .LT. -EPSO .OR. XI(2) .GT. TOL) GOTO 6100 + IF (XI(1)+XI(2) .LT. -EPSO .OR. XI(1)+XI(2) .GT. TOL) GOTO 6100 +C + GOTO 8000 +C +C Warning section +C + 6100 CONTINUE + IERR = 1 + IF (IPSW .GE. 0) THEN + WRITE(IWR,9961) XI,X + CALL DPRINT (TENC,NDIM,3,'TENC ',IWR) + ENDIF + GOTO 8000 +C +C Error section +C + 7200 CONTINUE + IERR = -2 + GOTO 7900 + 7300 CONTINUE + IERR = -3 + IF (IPSW .GE. 0) THEN + WRITE(IWR,9973) AREA + CALL DPRINT (TENC,NDIM,3,'TENC ',IWR) + ENDIF +C + 7900 CONTINUE + IF (IPSW .GE. 0) THEN + CALL ERRMSG ('T3_INV',IWR,IERR) + ENDIF +C +C +C Closing section +C + 8000 CONTINUE + IF (IPSW .GE. 7) WRITE(IWR,9910) + IF (IPSW .GE. 9) WRITE(IWR,9930) IERR,XI +C + RETURN +C + 9900 FORMAT( / ' ENTERING SUBROUTINE T3_INV' ) + 9910 FORMAT( / ' LEAVING SUBROUTINE T3_INV' ) + 9920 FORMAT( ' WITH INPUT PARAMETERS:' + + / ' NDIM = ', I6, + + / ' EPS = ', 1P,2E13.5, + + / ' X = ', 3E13.5 ) + 9930 FORMAT( ' WITH OUTPUT PARAMETERS:' + + / ' IERR = ', I6, + + / ' XI,ETA = ', 1P,2E13.5 ) + 9961 FORMAT( / ' ** WARNING FROM SUBROUTINE T3_INV' + + / ' THE POINT IS OUTSIDE THE ELEMENT', + + / ' XI,ETA = ', 1P,2E13.5, + + / ' X = ', 3E13.5 ) + 9973 FORMAT( / ' *** THE ELEMENT IS DEGENERATED', + + / ' AREA = ', 1P,E13.5 ) +C + END From fd550a8c399b443de346e11e12549afca76a6eab Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Thu, 6 Nov 2025 09:24:40 +0100 Subject: [PATCH 05/11] Added: virtual method ASMbase::findPoints() searching for a list of spatial points. The default implementation invokes ASMbase::findPoint() in a loop over the specified points. It is overridden in ASMu2DLag to perform parametric inversion for bilinear quads and triangles. --- src/ASM/ASMbase.C | 13 + src/ASM/ASMbase.h | 17 + src/ASM/ASMu2DLag.C | 193 +++ src/ASM/ASMu2DLag.h | 10 + src/ASM/Test/TestASMu2DLag.C | 73 + src/ASM/Test/refdata/Fine2Dmesh.xinp | 1861 ++++++++++++++++++++++++++ 6 files changed, 2167 insertions(+) create mode 100644 src/ASM/Test/refdata/Fine2Dmesh.xinp diff --git a/src/ASM/ASMbase.C b/src/ASM/ASMbase.C index 369f130b6..cd59c4563 100644 --- a/src/ASM/ASMbase.C +++ b/src/ASM/ASMbase.C @@ -1622,6 +1622,19 @@ bool ASMbase::injectNodeVec (const RealArray& nodeVec, RealArray& globRes, } +bool ASMbase::findPoints (std::vector& points, + std::vector& locs) const +{ + size_t nFound = 0; + locs.resize(points.size()); + for (size_t i = 0; i < points.size(); i++) + if ((locs[i].dist = this->findPoint(points[i],locs[i].u)) >= 0.0) + ++nFound; + + return nFound == points.size(); +} + + bool ASMbase::getSolution (Matrix& sField, const Vector& locSol, const IntVec& nodes) const { diff --git a/src/ASM/ASMbase.h b/src/ASM/ASMbase.h index 6aec0d545..9e0512aa7 100644 --- a/src/ASM/ASMbase.h +++ b/src/ASM/ASMbase.h @@ -582,6 +582,23 @@ class ASMbase //! \return Distance from the point \a X to the found point virtual double findPoint(Vec3& X, double* param) const = 0; + //! \brief A struct with element and associated local coordinates of a point. + struct PointParams + { + size_t iel; //!< 1-based element index + double u[3]; //!< Parametric coordinates w.r.t. element or patch + double dist; //!< Distance to spatial point + + //! \brief Default constructor. + explicit PointParams(size_t e = 0) : iel(e), u{0.0,0.0,0.0}, dist(-1.0) {} + }; + + //! \brief Finds elements and local coordinates matching spatial points. + //! \param points List of spatial point coordinates + //! \param[out] locs List of matching elements and/or parametric coordinates + virtual bool findPoints(std::vector& points, + std::vector& locs) const; + //! \brief Creates a standard FE model of this patch for visualization. //! \param[out] grid The generated finite element grid //! \param[in] npe Number of visualization nodes over each knot span diff --git a/src/ASM/ASMu2DLag.C b/src/ASM/ASMu2DLag.C index 5b96c2c0e..ec487651b 100644 --- a/src/ASM/ASMu2DLag.C +++ b/src/ASM/ASMu2DLag.C @@ -27,6 +27,24 @@ #include #endif +extern "C" { + //! \brief Parametric inversion for a linear triangle. + void t3_inv_(const int& nsd, const double& tol, const double& eps, + const double* X, const double* Xn, double* xi, + const int& ipsw, const int& iwr, int& ierr); + //! \brief Parametric inversion for a bilinear quadrilateral. + void q4_inv_(const int& nsd, const double& tol, const double& eps, + const double* X, const double* Xn, double* xi, + const int& ipsw, const int& iwr, int& ierr); + //! \brief Evaluate the bilinear shape functions. + void shape2_(const int& iop, const int& nenod, const double* xi, + double* shp, const int& ipsw, const int& iwr, int& ierr); + //! \brief Opens a temporary file for Fortran print. + int openftnfile_(const char* fname, const ssize_t nchar); + //! \brief Closes a Fortran file. + void closeftnfile_(const int& iunit); +} + ASMu2DLag::ASMu2DLag (unsigned char n_s, unsigned char n_f, char fType) : ASMs2DLag(n_s,n_f) @@ -619,3 +637,178 @@ bool ASMu2DLag::writeXML (const char* fname) const return true; } + + +bool ASMu2DLag::findPoints (Vec3Vec& points, ParamsVec& locs) const +{ + // Calculate the bounding box for all (shell) elements + size_t nonshell = 0; + std::vector bbox(nel); + for (size_t iel = 0; iel < nel; iel++) + if (MNPC[iel].size() == 3 || MNPC[iel].size() == 4) + this->getBoundingBox(MNPC[iel],bbox[iel]); + else + { + nonshell++; // make sure any other elements are sorted last + bbox[iel].first = bbox[iel].second = Vec3(1.0e99,1.0e99,1.0e99); + } + + std::array,3> eIdx; + std::vector::iterator it; + + // Index sort the elements in each spatial direction w.r.t. the max coordinate + for (int i = 0; i < 3; i++) + { + eIdx[i].resize(nel); + std::iota(eIdx[i].begin(), eIdx[i].end(), 0); + std::sort(eIdx[i].begin(), eIdx[i].end(), [i,&bbox](size_t a, size_t b) + { return bbox[a].second[i] < bbox[b].second[i]; }); + if (nonshell > 0) + eIdx[i].resize(nel-nonshell); // Remove the non-shell elements +#if SP_DEBUG > 2 + std::cout <<"\nElement order in "<< char('X'+i) <<"-direction:"; + for (size_t iel = 0; iel < nel-nonshell; iel++) + std::cout << (iel%10 ? ' ' : '\n') << MLGE[eIdx[i][iel]] + <<" ("<< bbox[eIdx[i][iel]].second[i] <<")"; + std::cout << std::endl; +#endif + } + + Vec3Pair BBox; // Bounding box for the whole patch + double maxDim = 0.0; + for (int i = 0; i < 3; i++) + { + BBox.first[i] = bbox[eIdx[i].front()].first[i]; + BBox.second[i] = bbox[eIdx[i].back()].second[i]; + if (double dim = BBox.second[i] - BBox.first[i]; dim > maxDim) + maxDim = dim; + } +#if SP_DEBUG > 1 + std::cout <<"\nDomain bounding box: [ " + << BBox.first <<" , "<< BBox.second + <<" ] maxDim = "<< maxDim << std::endl; +#endif + const double tol = maxDim*1.0e-8; + BBox.second += Vec3(RealArray(3,-tol)); + + size_t nFound = 0; + locs.clear(); + locs.reserve(points.size()); + for (Vec3& X : points) + { + // Do a binary search in each coordinate direction + // to find the range of elements to check in more detail + std::set nearElms; + for (int i = 0; i < 3; i++) + if (BBox.second[i] > BBox.first[i]) // skip the flat dimension + { + it = std::upper_bound(eIdx[i].begin(), eIdx[i].end(), X[i], + [i,&bbox](double x, size_t a) + { return x < bbox[a].second[i]; }); + if (nearElms.empty()) + nearElms.insert(it,eIdx[i].end()); + else + { + std::set tmp1(it,eIdx[i].end()), tmp2; + nearElms.swap(tmp2); + std::set_intersection(tmp1.begin(), tmp1.end(), + tmp2.begin(), tmp2.end(), + std::inserter(nearElms,nearElms.begin())); + } + } + + // Lambda function doing the Y < X comparison for two Vec3 objects + auto&& below = [&X,tol](const Vec3& Y) + { + for (int i = 0; i < 3; i++) + if (Y[i] > X[i]+tol) + return false; + return true; + }; + +#if SP_DEBUG > 1 + std::cout <<"\nSearching spatial point "<< X <<" (max "<< nearElms.size() + <<" elements to check)."<< std::endl; +#endif + bool found = false; + for (size_t iel : nearElms) + if (below(bbox[iel].first)) + { +#if SP_DEBUG > 1 + std::cout <<"Checking element "<< MLGE[iel] <<" in [ "<< bbox[iel].first + <<" - "<< bbox[iel].second <<" ]"<< std::endl; +#endif + // Perform parametric inversion to see if this point + // really is within current element + if (PointParams elem(1+iel); this->paramInvert(X,elem)) + { + locs.emplace_back(elem); + found = true; + break; + } + } + + if (found) + nFound++; + else + { + locs.emplace_back(PointParams()); + IFEM::cout <<" ** Warning: No element matches the point "<< X + << std::endl; + } + } + + return nFound == points.size(); +} + + +bool ASMu2DLag::paramInvert (Vec3& X, PointParams& elm) const +{ + Matrix Xn; + if (!this->getElementCoordinates(Xn,elm.iel)) + return false; + + const int nelnod = Xn.cols(); + if (nelnod == 4) // swap element node 3 and 4 + for (size_t r = 1; r <= Xn.rows(); r++) + std::swap(Xn(r,3),Xn(r,4)); + + const double eps = 1.0e-15; + const double tol = 1.0e-3; +#if SP_DEBUG > 2 + int ipsw = SP_DEBUG; +#else + int ipsw = 0; +#endif + // Open a temporary file for the Fortran output + std::string tmpfile = "/tmp/" + std::string(getenv("USER")) + ".ftn"; + const int iwr = openftnfile_(tmpfile.c_str(),tmpfile.size()); + if (iwr <= 0) ipsw = -1; // suppress all output + int ierr = 0; + if (nelnod == 3) + t3_inv_(nsd,tol,eps,X.ptr(),Xn.ptr(),elm.u,ipsw,iwr,ierr); + else if (nelnod == 4) + q4_inv_(nsd,tol,eps,X.ptr(),Xn.ptr(),elm.u,ipsw,iwr,ierr); + else + ierr = -99; + if (iwr > 0 && ipsw >= 0) + { + // Copy fortran output to IFEM::cout + closeftnfile_(iwr); + std::string cline; + std::ifstream is(tmpfile); + while (std::getline(is,cline)) + IFEM::cout << cline << std::endl; + } + if (ierr != 0 && ierr != 2) + return false; + + // Calculate the distance to the projected point + Vector N(nelnod); + shape2_(0,nelnod,elm.u,N.ptr(),ipsw,iwr,ierr); + Vec3 X0(Xn * N); + elm.dist = (X0-X).length(); + + X = X0; // update the spatial coordinates + return true; +} diff --git a/src/ASM/ASMu2DLag.h b/src/ASM/ASMu2DLag.h index 6b30ced43..3fb75abb2 100644 --- a/src/ASM/ASMu2DLag.h +++ b/src/ASM/ASMu2DLag.h @@ -136,6 +136,13 @@ 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; + using ParamsVec = std::vector; //!< Convenience type alias + + //! \brief Finds elements and local coordinates matching spatial points. + //! \param points List of spatial point coordinates + //! \param[out] locs List of matching element and local coordinates + virtual bool findPoints(Vec3Vec& points, ParamsVec& locs) const; + //! \brief Dumps the mesh to the specified XML-file. bool writeXML(const char* fname) const; @@ -143,6 +150,9 @@ class ASMu2DLag : public ASMs2DLag //! \brief Generate thread groups using multi-coloring. void generateThreadGroupsMultiColored(bool silence, bool separateGroup1Noded); + //! \brief Performs parametric inversion for an element. + bool paramInvert(Vec3& X, PointParams& elm) const; + bool swapNode34; //!< If \e true, element nodes 3 and 4 should be swapped std::vector nodeSets; //!< Node sets for Dirichlet BCs diff --git a/src/ASM/Test/TestASMu2DLag.C b/src/ASM/Test/TestASMu2DLag.C index c045e8ae9..136e67e2d 100644 --- a/src/ASM/Test/TestASMu2DLag.C +++ b/src/ASM/Test/TestASMu2DLag.C @@ -13,7 +13,9 @@ #include "Catch2Support.h" #include "ASMu2DLag.h" #include "MPC.h" +#include "Vec3Oper.h" +#include #include @@ -224,3 +226,74 @@ TEST_CASE("TestASMu2DLag.GenerateThreadGroups3x3OneNodeSPC") } } } + + +TEST_CASE("TestASMu2DLag.findPoints") +{ + struct TestCase + { + const char* fname = nullptr; + int nX = 0; + int nY = 0; + size_t npt = 0; + Vec3 X; + }; + + const TestCase param = + GENERATE(TestCase{ "5x4 mesh", 5, 4, 3, Vec3(0.123,0.345,0.0) }, + TestCase{ "src/ASM/Test/refdata/Fine2Dmesh.xinp", 0, 0, 17, + Vec3(0.1,0.0,0.0) } + ); + + SECTION(param.fname) + { + std::cout <<"\nTesting "<< param.fname; + if (param.nX > 0) std::cout <<" "<< param.nX <<"x"<< param.nY; + std::cout <<" "<< param.npt << std::endl; + + ASMbase::resetNumbering(); + ASMu2DLag pch(3,3,'x'); + if (param.nX > 0 && param.nY > 0) + { + // Generate a simple mesh on the bi-unit square + std::stringstream str; + REQUIRE(pch.read(generateXMLModel(str,1.0,1.0,param.nX,param.nY))); + } + else + { + // Read mesh from file + std::ifstream fs(param.fname); + REQUIRE(pch.read(fs)); + } + REQUIRE(pch.generateFEMTopology()); + + // Pick some random elements and search for their center points + const size_t nel = pch.getNoElms(); + const size_t npt = param.npt; + std::vector elements; elements.reserve(npt); + std::vector points; points.reserve(npt); + for (size_t iel = 1+nel/(2*npt); iel <= nel; iel += nel/npt) + { + elements.push_back(iel); + points.emplace_back(pch.getElementCenter(iel)); + std::cout <<"Element "<< iel <<": "<< points.back() << std::endl; + } + points.emplace_back(param.X); + + std::vector elms; + REQUIRE(pch.findPoints(points,elms)); + + size_t ipt = 0; + std::cout <<"Search results:"; + for (const ASMbase::PointParams& elm : elms) + { + std::cout <<"\nPoint "<< ++ipt <<": iel = "<< elm.iel + <<" xi,eta = "<< elm.u[0] <<","<< elm.u[1] + <<" dist = "<< elm.dist; + if (ipt-1 < elements.size()) + REQUIRE(elements[ipt-1] == elm.iel); + REQUIRE(elm.dist <= 1.0e-12); + } + std::cout << std::endl; + } +} diff --git a/src/ASM/Test/refdata/Fine2Dmesh.xinp b/src/ASM/Test/refdata/Fine2Dmesh.xinp new file mode 100644 index 000000000..8f5c832b2 --- /dev/null +++ b/src/ASM/Test/refdata/Fine2Dmesh.xinp @@ -0,0 +1,1861 @@ + + + 0.100488 0 -0.00270879 + 0.100458 0 -0.00246075 + 0.100423 0 -0.00226954 + 0.100375 0 -0.00207408 + 0.100287 0 -0.00179642 + 0.100175 0 -0.00152714 + 0.100041 0 -0.00126816 + 0.0999368 0 -0.00109605 + 0.0998266 0 -0.000936094 + 0.0996721 0 -0.00073968 + 0.0994749 0 -0.000525124 + 0.0992603 0 -0.000327851 + 0.0990639 0 -0.000173401 + 0.0989039 0 -6.31507e-05 + 0.0987318 0 4.12858e-05 + 0.0984729 0 0.000174987 + 0.0982036 0 0.000286524 + 0.0979259 0 0.000375111 + 0.0977304 0 0.000422931 + 0.0975393 0 0.000458118 + 0.0972912 0 0.000487782 + 0.100987 0 -0.00267906 + 0.100954 0 -0.00239547 + 0.100911 0 -0.00216132 + 0.100852 0 -0.00192463 + 0.100754 0 -0.00161889 + 0.100631 0 -0.00132213 + 0.100484 0 -0.00103632 + 0.100359 0 -0.000827405 + 0.100223 0 -0.00063166 + 0.100046 0 -0.000407722 + 0.0998284 0 -0.000171571 + 0.0995923 0 4.61626e-05 + 0.0993683 0 0.00022334 + 0.0991726 0 0.000358553 + 0.0989637 0 0.000484495 + 0.0986779 0 0.000631049 + 0.0983811 0 0.00075397 + 0.0980754 0 0.000852446 + 0.0978386 0 0.000911088 + 0.0976046 0 0.000953922 + 0.0973209 0 0.000986979 + 0.101486 0 -0.00264932 + 0.10145 0 -0.00233019 + 0.101399 0 -0.00205311 + 0.10133 0 -0.00177517 + 0.101221 0 -0.00144137 + 0.101087 0 -0.00111712 + 0.100928 0 -0.000804473 + 0.10078 0 -0.000558757 + 0.10062 0 -0.000327226 + 0.10042 0 -7.57639e-05 + 0.100182 0 0.000181982 + 0.0999242 0 0.000420176 + 0.0996728 0 0.000620082 + 0.0994412 0 0.000780257 + 0.0991955 0 0.000927705 + 0.0988829 0 0.00108711 + 0.0985586 0 0.00122142 + 0.0982248 0 0.00132978 + 0.0979468 0 0.00139924 + 0.0976698 0 0.00144973 + 0.0973507 0 0.00148618 + 0.101985 0 -0.00261958 + 0.101946 0 -0.00226491 + 0.101887 0 -0.0019449 + 0.101807 0 -0.00162571 + 0.101689 0 -0.00126384 + 0.101543 0 -0.000912119 + 0.101371 0 -0.00057263 + 0.101202 0 -0.000290108 + 0.101017 0 -2.27914e-05 + 0.100794 0 0.000256194 + 0.100536 0 0.000535535 + 0.100256 0 0.00079419 + 0.0999772 0 0.00101682 + 0.0997099 0 0.00120196 + 0.0994274 0 0.00137091 + 0.0990879 0 0.00154317 + 0.0987362 0 0.00168886 + 0.0983743 0 0.00180711 + 0.0980551 0 0.0018874 + 0.0977351 0 0.00194553 + 0.0973804 0 0.00198537 + 0.102485 0 -0.00258984 + 0.102441 0 -0.00219963 + 0.102376 0 -0.00183668 + 0.102284 0 -0.00147625 + 0.102156 0 -0.00108632 + 0.101999 0 -0.000707113 + 0.101814 0 -0.000340786 + 0.101624 0 -2.14593e-05 + 0.101414 0 0.000281643 + 0.101168 0 0.000588152 + 0.100889 0 0.000889088 + 0.100588 0 0.0011682 + 0.100282 0 0.00141356 + 0.0999785 0 0.00162366 + 0.0996592 0 0.00181412 + 0.0992929 0 0.00199924 + 0.0989137 0 0.00215631 + 0.0985238 0 0.00228445 + 0.0981633 0 0.00237556 + 0.0978004 0 0.00244134 + 0.0974102 0 0.00248457 + 0.102984 0 -0.0025601 + 0.102937 0 -0.00213435 + 0.102864 0 -0.00172847 + 0.102762 0 -0.00132679 + 0.102624 0 -0.000908797 + 0.102455 0 -0.000502107 + 0.102257 0 -0.000108943 + 0.102045 0 0.000247189 + 0.10181 0 0.000586077 + 0.101542 0 0.00092011 + 0.101243 0 0.00124264 + 0.10092 0 0.00154222 + 0.100586 0 0.00181031 + 0.100247 0 0.00204537 + 0.0998911 0 0.00225733 + 0.0994979 0 0.0024553 + 0.0990912 0 0.00262376 + 0.0986732 0 0.00276178 + 0.0982715 0 0.00286371 + 0.0978657 0 0.00293714 + 0.0974399 0 0.00298377 + 0.103483 0 -0.00253036 + 0.103433 0 -0.00206908 + 0.103352 0 -0.00162026 + 0.103239 0 -0.00117733 + 0.103091 0 -0.000731274 + 0.102911 0 -0.000297101 + 0.102701 0 0.000122901 + 0.102467 0 0.000515838 + 0.102207 0 0.000890511 + 0.101916 0 0.00125207 + 0.101596 0 0.00159619 + 0.101252 0 0.00191623 + 0.100891 0 0.00220705 + 0.100516 0 0.00246707 + 0.100123 0 0.00270054 + 0.0997029 0 0.00291136 + 0.0992687 0 0.0030912 + 0.0988227 0 0.00323912 + 0.0983797 0 0.00335187 + 0.0979309 0 0.00343295 + 0.0974696 0 0.00348297 + 0.097 0 0.0035 + 0.097 0 0.003 + 0.097 0 0.0025 + 0.097 0 0.002 + 0.097 0 0.0015 + 0.097 0 0.001 + 0.097 0 0.0005 + 0.0972615 0 -1.14159e-05 + 0.097474 0 -3.76868e-05 + 0.0976222 0 -6.52253e-05 + 0.0977765 0 -0.000102223 + 0.0980261 0 -0.000180923 + 0.0982679 0 -0.000281076 + 0.0985 0 -0.000401924 + 0.0986352 0 -0.000484854 + 0.0987595 0 -0.000570142 + 0.0989284 0 -0.000701865 + 0.0991213 0 -0.000878677 + 0.0992981 0 -0.00107164 + 0.0994299 0 -0.00124053 + 0.0995151 0 -0.0013647 + 0.0995981 0 -0.0015 + 0.0997189 0 -0.00173214 + 0.0998191 0 -0.00197394 + 0.0998978 0 -0.00222354 + 0.0999348 0 -0.00237775 + 0.0999623 0 -0.00252603 + 0.0999886 0 -0.00273853 + 0.103982 0 -0.00250063 + 0.103929 0 -0.0020038 + 0.10384 0 -0.00151204 + 0.103716 0 -0.00102787 + 0.103559 0 -0.00055375 + 0.103367 0 -9.20946e-05 + 0.103144 0 0.000354744 + 0.102889 0 0.000784487 + 0.102604 0 0.00119495 + 0.10229 0 0.00158403 + 0.10195 0 0.00194975 + 0.101584 0 0.00229025 + 0.101195 0 0.00260379 + 0.100784 0 0.00288877 + 0.100355 0 0.00314375 + 0.0999079 0 0.00336742 + 0.0994462 0 0.00355865 + 0.0989721 0 0.00371645 + 0.098488 0 0.00384003 + 0.0979962 0 0.00392875 + 0.0974994 0 0.00398216 + 0.1035 0 -0.003 + 0.103 0 -0.003 + 0.1025 0 -0.003 + 0.102 0 -0.003 + 0.1015 0 -0.003 + 0.101 0 -0.003 + 0.1005 0 -0.003 + 0.097 0 0 + 0.097 0 0.004 + 0.1 0 -0.003 + 0.104 0 -0.003 + 0.12665 0 0.0289971 + 0.125661 0 0.0289445 + 0.124446 0 0.0288412 + 0.123023 0 0.0287277 + 0.12148 0 0.0285904 + 0.119865 0 0.0283934 + 0.11806 0 0.0283224 + 0.116262 0 0.0282936 + 0.114485 0 0.0282789 + 0.112723 0 0.0282689 + 0.110969 0 0.0282602 + 0.10922 0 0.0282521 + 0.107472 0 0.0282442 + 0.105726 0 0.0282367 + 0.10398 0 0.0282296 + 0.102235 0 0.0282231 + 0.10049 0 0.0282171 + 0.0987449 0 0.0282118 + 0.0987405 0 0.026423 + 0.0987342 0 0.0246333 + 0.0987252 0 0.0228425 + 0.0987135 0 0.0210504 + 0.0986995 0 0.0192565 + 0.0986837 0 0.0174593 + 0.0986664 0 0.0156552 + 0.0986454 0 0.0138387 + 0.0986115 0 0.0120046 + 0.09853 0 0.0101753 + 0.098273 0 0.00856521 + 0.0980445 0 0.00722559 + 0.0978595 0 0.0061374 + 0.0977098 0 0.00527673 + 0.0975947 0 0.00457822 + 0.0981704 0 0.00452173 + 0.0987063 0 0.00439854 + 0.0991867 0 0.00425457 + 0.0996818 0 0.00406872 + 0.100156 0 0.00385058 + 0.100595 0 0.003631 + 0.101063 0 0.00337721 + 0.101526 0 0.00307077 + 0.101952 0 0.00272058 + 0.102352 0 0.00235359 + 0.102687 0 0.00194511 + 0.102997 0 0.00155958 + 0.103329 0 0.00114126 + 0.10363 0 0.000666665 + 0.103862 0 0.000173918 + 0.104053 0 -0.000317079 + 0.104215 0 -0.000798186 + 0.104378 0 -0.00128132 + 0.104517 0 -0.00181102 + 0.104579 0 -0.0023927 + 0.105284 0 -0.00225575 + 0.106168 0 -0.0020612 + 0.107303 0 -0.00180852 + 0.108726 0 -0.00169339 + 0.110311 0 -0.00155966 + 0.11205 0 -0.00146469 + 0.113841 0 -0.00140362 + 0.115644 0 -0.00136419 + 0.117449 0 -0.00133748 + 0.11925 0 -0.0013181 + 0.121049 0 -0.00130295 + 0.122844 0 -0.00129034 + 0.124636 0 -0.00127939 + 0.126426 0 -0.00126962 + 0.128213 0 -0.00126058 + 0.12822 0 0.00047981 + 0.128226 0 0.00222217 + 0.128231 0 0.00396681 + 0.128236 0 0.00571375 + 0.128241 0 0.00746276 + 0.128246 0 0.0092137 + 0.128252 0 0.010967 + 0.12826 0 0.0127245 + 0.12827 0 0.0144899 + 0.128286 0 0.01627 + 0.128319 0 0.0180724 + 0.128397 0 0.0198807 + 0.128602 0 0.0214944 + 0.128746 0 0.023027 + 0.128852 0 0.0244135 + 0.128902 0 0.025545 + 0.12891 0 0.0265387 + 0.128815 0 0.0274124 + 0.128556 0 0.0281426 + 0.128133 0 0.0286724 + 0.127481 0 0.0289285 + 0.126466 0 0.0281245 + 0.12554 0 0.0279441 + 0.124385 0 0.027695 + 0.123013 0 0.0274789 + 0.121569 0 0.0272373 + 0.120205 0 0.026656 + 0.118144 0 0.0265956 + 0.116277 0 0.0265654 + 0.114476 0 0.0265459 + 0.112705 0 0.0265298 + 0.110947 0 0.0265145 + 0.109197 0 0.0264994 + 0.10745 0 0.0264846 + 0.105706 0 0.0264703 + 0.103963 0 0.0264568 + 0.102222 0 0.0264444 + 0.100481 0 0.0264331 + 0.100469 0 0.0246475 + 0.100452 0 0.0228596 + 0.10043 0 0.0210693 + 0.100405 0 0.0192759 + 0.100377 0 0.0174767 + 0.100348 0 0.0156658 + 0.100318 0 0.0138297 + 0.100285 0 0.0119321 + 0.100249 0 0.00985007 + 0.0995275 0 0.00824778 + 0.0990525 0 0.00698662 + 0.0986878 0 0.00595154 + 0.0983874 0 0.00518214 + 0.0989839 0 0.00497915 + 0.0993911 0 0.00483662 + 0.0999436 0 0.00461402 + 0.100445 0 0.0043389 + 0.100815 0 0.00415674 + 0.101353 0 0.00392321 + 0.101903 0 0.00358705 + 0.102354 0 0.00317382 + 0.102826 0 0.00280539 + 0.103114 0 0.00228998 + 0.103373 0 0.00196416 + 0.103806 0 0.00156255 + 0.104187 0 0.00100622 + 0.1044 0 0.000448633 + 0.10458 0 -7.92453e-05 + 0.104715 0 -0.000555746 + 0.104944 0 -0.000994438 + 0.105183 0 -0.00155905 + 0.10599 0 -0.00116758 + 0.10711 0 -0.000464572 + 0.108683 0 -0.000390918 + 0.110195 0 -6.73027e-05 + 0.111985 0 0.000115755 + 0.113812 0 0.000223629 + 0.115636 0 0.000292083 + 0.117452 0 0.000339106 + 0.119261 0 0.000374 + 0.121063 0 0.000401768 + 0.122859 0 0.000425108 + 0.12465 0 0.000445418 + 0.126436 0 0.000463366 + 0.126447 0 0.00219949 + 0.126457 0 0.00393984 + 0.126467 0 0.00568469 + 0.126476 0 0.00743383 + 0.126487 0 0.00918698 + 0.126499 0 0.0109447 + 0.126512 0 0.0127094 + 0.126529 0 0.0144874 + 0.126552 0 0.0162936 + 0.126591 0 0.0181674 + 0.12667 0 0.0202362 + 0.127266 0 0.0216007 + 0.127533 0 0.0230384 + 0.127761 0 0.0243603 + 0.127849 0 0.0252196 + 0.127923 0 0.0262034 + 0.127939 0 0.0270476 + 0.12785 0 0.0277287 + 0.127732 0 0.0281975 + 0.127213 0 0.0281954 + 0.123089 0 0.0262564 + 0.124556 0 0.0265196 + 0.127174 0 0.0275357 + 0.102204 0 0.0246633 + 0.103941 0 0.0246809 + 0.105681 0 0.0247 + 0.107423 0 0.0247201 + 0.109168 0 0.0247408 + 0.110917 0 0.0247613 + 0.112674 0 0.024781 + 0.114444 0 0.0247995 + 0.116233 0 0.0248165 + 0.118048 0 0.02483 + 0.119872 0 0.0248264 + 0.121525 0 0.0247355 + 0.123211 0 0.0248122 + 0.125098 0 0.025117 + 0.102181 0 0.022879 + 0.103914 0 0.0229008 + 0.10565 0 0.0229245 + 0.10739 0 0.0229495 + 0.109133 0 0.0229749 + 0.110882 0 0.0229998 + 0.112637 0 0.0230232 + 0.114401 0 0.0230443 + 0.116175 0 0.0230624 + 0.117956 0 0.0230766 + 0.119726 0 0.0230846 + 0.121441 0 0.0230872 + 0.123143 0 0.0231488 + 0.124838 0 0.0232475 + 0.102152 0 0.0210915 + 0.10388 0 0.0211168 + 0.105613 0 0.0211445 + 0.107351 0 0.0211736 + 0.109094 0 0.0212029 + 0.110842 0 0.0212316 + 0.112596 0 0.0212586 + 0.114356 0 0.0212833 + 0.116122 0 0.0213056 + 0.117888 0 0.021327 + 0.119645 0 0.0213511 + 0.121378 0 0.0213876 + 0.123089 0 0.0214589 + 0.124756 0 0.0215609 + 0.102119 0 0.0193003 + 0.103841 0 0.019329 + 0.105571 0 0.0193606 + 0.107308 0 0.0193935 + 0.10905 0 0.0194266 + 0.110799 0 0.0194588 + 0.112552 0 0.0194894 + 0.114311 0 0.0195182 + 0.116074 0 0.0195457 + 0.117837 0 0.0195745 + 0.119595 0 0.0196099 + 0.121342 0 0.0196621 + 0.123083 0 0.0197495 + 0.124833 0 0.0199048 + 0.102082 0 0.0175034 + 0.103799 0 0.0175367 + 0.105525 0 0.0175731 + 0.107261 0 0.0176104 + 0.109003 0 0.0176472 + 0.110752 0 0.0176829 + 0.112507 0 0.0177171 + 0.114267 0 0.01775 + 0.11603 0 0.0177824 + 0.117794 0 0.0178167 + 0.119556 0 0.0178571 + 0.121313 0 0.0179098 + 0.123066 0 0.0179822 + 0.124823 0 0.0180772 + 0.102045 0 0.0156972 + 0.103755 0 0.0157387 + 0.105478 0 0.0157823 + 0.107211 0 0.015825 + 0.108953 0 0.0158659 + 0.110702 0 0.0159051 + 0.112459 0 0.0159429 + 0.114221 0 0.0159796 + 0.115987 0 0.0160162 + 0.117754 0 0.0160542 + 0.119521 0 0.0160961 + 0.121284 0 0.0161445 + 0.123043 0 0.0162001 + 0.124799 0 0.0162578 + 0.10201 0 0.0138763 + 0.103713 0 0.0139345 + 0.105429 0 0.0139892 + 0.107158 0 0.0140387 + 0.108899 0 0.014084 + 0.110649 0 0.0141265 + 0.112408 0 0.0141675 + 0.114173 0 0.0142077 + 0.115943 0 0.0142479 + 0.117715 0 0.0142888 + 0.119486 0 0.0143314 + 0.121255 0 0.0143759 + 0.123018 0 0.0144212 + 0.124776 0 0.0144624 + 0.101982 0 0.012038 + 0.103675 0 0.0121291 + 0.10538 0 0.0121979 + 0.107103 0 0.0122539 + 0.108841 0 0.0123028 + 0.110591 0 0.012348 + 0.112353 0 0.0123917 + 0.114122 0 0.0124351 + 0.115897 0 0.0124786 + 0.117674 0 0.0125221 + 0.119451 0 0.0125654 + 0.121225 0 0.0126081 + 0.122994 0 0.0126486 + 0.124756 0 0.0126842 + 0.101979 0 0.0102086 + 0.10364 0 0.0103413 + 0.105327 0 0.0104162 + 0.107042 0 0.0104743 + 0.108776 0 0.0105248 + 0.110528 0 0.0105711 + 0.112293 0 0.0106168 + 0.114068 0 0.010663 + 0.115849 0 0.0107096 + 0.117633 0 0.0107556 + 0.119417 0 0.0108002 + 0.121197 0 0.0108426 + 0.122971 0 0.010882 + 0.124737 0 0.0109167 + 0.102064 0 0.00860037 + 0.103598 0 0.00860859 + 0.105262 0 0.00864992 + 0.106971 0 0.00870297 + 0.108703 0 0.0087527 + 0.110455 0 0.00879785 + 0.112226 0 0.00884427 + 0.114009 0 0.00889317 + 0.1158 0 0.00894289 + 0.117592 0 0.00899148 + 0.119384 0 0.0090377 + 0.12117 0 0.00908094 + 0.122949 0 0.00912069 + 0.124721 0 0.00915629 + 0.101756 0 0.0069012 + 0.103443 0 0.00684337 + 0.105163 0 0.00687461 + 0.106886 0 0.00693964 + 0.108616 0 0.00699155 + 0.110369 0 0.00703038 + 0.11215 0 0.00707616 + 0.113947 0 0.0071283 + 0.11575 0 0.00718182 + 0.117553 0 0.00723325 + 0.119352 0 0.00728114 + 0.121145 0 0.0073252 + 0.122929 0 0.00736544 + 0.124705 0 0.00740179 + 0.101965 0 0.00546758 + 0.103269 0 0.0049944 + 0.105072 0 0.00507335 + 0.106804 0 0.00519968 + 0.108513 0 0.00525533 + 0.11026 0 0.00526811 + 0.112064 0 0.00531278 + 0.113883 0 0.00537142 + 0.115702 0 0.0054306 + 0.117516 0 0.00548516 + 0.119322 0 0.00553414 + 0.12112 0 0.00557806 + 0.12291 0 0.00561762 + 0.12469 0 0.00565332 + 0.103493 0 0.00341329 + 0.105061 0 0.00323793 + 0.106754 0 0.00354626 + 0.108377 0 0.00357906 + 0.110101 0 0.00349043 + 0.111968 0 0.00354996 + 0.113825 0 0.00362628 + 0.115662 0 0.00369438 + 0.117485 0 0.00375168 + 0.119296 0 0.00380007 + 0.121098 0 0.00384182 + 0.122891 0 0.00387858 + 0.124676 0 0.00391144 + 0.105658 0 0.00203031 + 0.106778 0 0.00218835 + 0.108146 0 0.00204454 + 0.109806 0 0.0015831 + 0.111887 0 0.00178606 + 0.11379 0 0.00190254 + 0.115638 0 0.00198019 + 0.117463 0 0.0020369 + 0.119275 0 0.00208118 + 0.121078 0 0.00211761 + 0.122873 0 0.00214883 + 0.124661 0 0.00217629 + 0.101643 0 0.00457787 + 0.100793 0 0.00600991 + 0.100265 0 0.00521521 + 0.101209 0 0.00540471 + 0.100951 0 0.00867966 + 0.09996 0 0.00652133 + 0.100549 0 0.00758749 + 0.102744 0 0.00359039 + 0.10434 0 0.0021509 + 0.107629 0 0.000848334 + 0.105156 0 0.000123113 + 0.10595 0 0.000225313 + 0.106559 0 0.00115345 + 0.126361 0 0.0231743 + 0.10239 0 0.00418743 + 0.105882 0 0.00136163 + 0.105124 0 -0.000336288 + 0.105502 0 -0.000568249 + 0.126466 0 0.0273734 + 0.12566 0 0.0270188 + 0.108706 0 0.000680171 + 0.126995 0 0.0260161 + 0.127169 0 0.0268546 + 0.126081 0 0.0262427 + 0.126572 0 0.0268236 + 0.121597 0 0.0262213 + 0.100875 0 0.00473868 + 0.0994576 0 0.00550077 + 0.10358 0 0.00245458 + 0.10543 0 0.000910276 + 0.104975 0 0.000708883 + 0.104918 0 0.00136132 + 0.126263 0 0.021643 + 0.126813 0 0.0247775 + 0.097 0 0.0282067 + 0.097 0 0.0264134 + 0.097 0 0.0246201 + 0.097 0 0.0228269 + 0.097 0 0.0210336 + 0.097 0 0.0192403 + 0.097 0 0.017447 + 0.097 0 0.0156537 + 0.097 0 0.0138604 + 0.097 0 0.0120672 + 0.097 0 0.0102772 + 0.097 0 0.00860989 + 0.097 0 0.00721229 + 0.097 0 0.00609581 + 0.097 0 0.00520954 + 0.097 0 0.0045325 + 0.126014 0 0.03 + 0.124727 0 0.03 + 0.123168 0 0.03 + 0.121476 0 0.03 + 0.119728 0 0.03 + 0.117979 0 0.03 + 0.116231 0 0.03 + 0.114483 0 0.03 + 0.112735 0 0.03 + 0.110986 0 0.03 + 0.109238 0 0.03 + 0.10749 0 0.03 + 0.105741 0 0.03 + 0.103993 0 0.03 + 0.102245 0 0.03 + 0.100497 0 0.03 + 0.0987483 0 0.03 + 0.128207 0 -0.003 + 0.126413 0 -0.003 + 0.12462 0 -0.003 + 0.122827 0 -0.003 + 0.121034 0 -0.003 + 0.11924 0 -0.003 + 0.117447 0 -0.003 + 0.115654 0 -0.003 + 0.11386 0 -0.003 + 0.112067 0 -0.003 + 0.110277 0 -0.003 + 0.10861 0 -0.003 + 0.107212 0 -0.003 + 0.106096 0 -0.003 + 0.10521 0 -0.003 + 0.104533 0 -0.003 + 0.129853 0 0.027927 + 0.129427 0 0.0287634 + 0.128763 0 0.029427 + 0.127927 0 0.0298532 + 0.13 0 -0.00125171 + 0.13 0 0.000496575 + 0.13 0 0.00224486 + 0.13 0 0.00399315 + 0.13 0 0.00574144 + 0.13 0 0.00748972 + 0.13 0 0.00923801 + 0.13 0 0.0109863 + 0.13 0 0.0127346 + 0.13 0 0.0144829 + 0.13 0 0.0162312 + 0.13 0 0.0179794 + 0.13 0 0.0197277 + 0.13 0 0.021476 + 0.13 0 0.0231677 + 0.13 0 0.0247271 + 0.13 0 0.0260137 + 0.097 0 0.03 + 0.127 0 0.03 + 0.13 0 -0.003 + 0.13 0 0.027 + 0.0884378 0 0.0281585 + 0.0901539 0 0.0281953 + 0.0918678 0 0.0282051 + 0.09358 0 0.0282077 + 0.0952908 0 0.0282077 + 0.0884483 0 0.0264379 + 0.0901657 0 0.0264177 + 0.0918785 0 0.0264174 + 0.0935878 0 0.0264178 + 0.0952946 0 0.0264163 + 0.0867357 0 0.0245771 + 0.0884638 0 0.024603 + 0.0901823 0 0.0246203 + 0.0918929 0 0.0246288 + 0.093598 0 0.02463 + 0.0952998 0 0.0246263 + 0.0884889 0 0.0227769 + 0.090207 0 0.0228318 + 0.0919128 0 0.0228475 + 0.0936115 0 0.0228471 + 0.0953065 0 0.0228388 + 0.088532 0 0.021085 + 0.0902441 0 0.0210823 + 0.0919399 0 0.0210821 + 0.0936287 0 0.0210722 + 0.0953148 0 0.0210547 + 0.0868662 0 0.0192233 + 0.0886073 0 0.0193021 + 0.0902975 0 0.0193303 + 0.091974 0 0.0193264 + 0.0936487 0 0.0193049 + 0.0953241 0 0.0192744 + 0.0870697 0 0.0173217 + 0.0887337 0 0.0175696 + 0.0903646 0 0.0176103 + 0.09201 0 0.0175883 + 0.0936681 0 0.0175464 + 0.0953328 0 0.0174977 + 0.0876791 0 0.0161301 + 0.0888933 0 0.0160442 + 0.0904174 0 0.0159531 + 0.0920335 0 0.0158701 + 0.0936809 0 0.0157944 + 0.095339 0 0.0157232 + 0.0870745 0 0.0150243 + 0.0887432 0 0.0145238 + 0.0903782 0 0.0142874 + 0.0920259 0 0.0141444 + 0.0936831 0 0.0140381 + 0.0953424 0 0.0139468 + 0.086876 0 0.0130795 + 0.0886271 0 0.012739 + 0.0903264 0 0.012528 + 0.092009 0 0.0123817 + 0.0936832 0 0.0122666 + 0.0953472 0 0.012166 + 0.0868032 0 0.0109182 + 0.0885635 0 0.0108242 + 0.0902921 0 0.0107037 + 0.0920009 0 0.0105875 + 0.0936922 0 0.0104828 + 0.0953618 0 0.0103857 + 0.0885347 0 0.00893657 + 0.0902799 0 0.00887591 + 0.0920082 0 0.00878814 + 0.0937173 0 0.00870094 + 0.0954037 0 0.0086223 + 0.0885292 0 0.0071904 + 0.090286 0 0.0070862 + 0.0920265 0 0.00700692 + 0.0937526 0 0.00692854 + 0.0955286 0 0.0068007 + 0.0867698 0 0.00540108 + 0.0885411 0 0.00534604 + 0.0903062 0 0.00529404 + 0.0920444 0 0.00525206 + 0.0937195 0 0.0052281 + 0.0952284 0 0.00522512 + 0.0885675 0 0.00351718 + 0.090348 0 0.00351661 + 0.0921082 0 0.00350542 + 0.0938315 0 0.00352788 + 0.0954514 0 0.00369025 + 0.0886041 0 0.00181618 + 0.0904041 0 0.00176717 + 0.0921952 0 0.00174188 + 0.0960787 0 0.00147056 + 0.0965341 0 0.00184952 + 0.0963531 0 0.00243252 + 0.0961985 0 0.00519451 + 0.0963132 0 0.00432189 + 0.0962424 0 0.00603034 + 0.0964519 0 0.000612568 + 0.095818 0 0.000807326 + 0.0965991 0 0.00387286 + 0.0965576 0 0.00114515 + 0.0966972 0 0.00149822 + 0.0867196 0 0.0280008 + 0.0867259 0 0.0265721 + 0.0867531 0 0.0225877 + 0.0867878 0 0.0211785 + 0.0867758 0 0.00886077 + 0.0867687 0 0.00740852 + 0.0867846 0 0.00340786 + 0.0868026 0 0.00198783 + 0.0940323 0 0.00169997 + 0.0953878 0 0.0020879 + 0.0950197 0 0.0011488 + 0.0964807 0 0.00329913 + 0.103638 0 -0.0046441 + 0.105189 0 -0.00481537 + 0.106779 0 -0.00449506 + 0.108593 0 -0.00461211 + 0.11034 0 -0.0046468 + 0.112103 0 -0.00465749 + 0.113872 0 -0.00466331 + 0.115648 0 -0.00467149 + 0.117432 0 -0.00468258 + 0.119225 0 -0.00469506 + 0.121024 0 -0.00470759 + 0.122825 0 -0.00471919 + 0.124627 0 -0.0047284 + 0.126426 0 -0.00473243 + 0.128215 0 -0.00472564 + 0.103509 0 -0.00622966 + 0.105205 0 -0.00633575 + 0.106897 0 -0.00628743 + 0.108644 0 -0.0063066 + 0.110387 0 -0.00631763 + 0.112132 0 -0.00631988 + 0.113877 0 -0.00632435 + 0.115633 0 -0.00634014 + 0.117409 0 -0.00636385 + 0.119204 0 -0.00639016 + 0.121011 0 -0.00641622 + 0.122823 0 -0.00644088 + 0.124637 0 -0.00646212 + 0.126439 0 -0.00646766 + 0.128223 0 -0.00645166 + 0.101737 0 -0.00783441 + 0.103493 0 -0.00794336 + 0.105231 0 -0.00800666 + 0.106963 0 -0.00801005 + 0.108699 0 -0.00800812 + 0.110434 0 -0.0079968 + 0.11216 0 -0.00797995 + 0.113872 0 -0.00797409 + 0.115598 0 -0.00800096 + 0.117367 0 -0.00804261 + 0.119172 0 -0.00808571 + 0.120993 0 -0.0081267 + 0.122818 0 -0.00816657 + 0.124644 0 -0.0081964 + 0.126447 0 -0.00819828 + 0.128224 0 -0.00817363 + 0.10176 0 -0.00962154 + 0.103504 0 -0.00969011 + 0.105269 0 -0.00972834 + 0.107028 0 -0.00973266 + 0.108756 0 -0.00971642 + 0.110493 0 -0.00968036 + 0.112202 0 -0.00962863 + 0.113854 0 -0.00959102 + 0.115521 0 -0.00964705 + 0.117288 0 -0.00972021 + 0.119125 0 -0.00978428 + 0.120974 0 -0.00983971 + 0.1228 0 -0.00988831 + 0.124636 0 -0.00991727 + 0.126447 0 -0.00991438 + 0.128213 0 -0.00988758 + 0.101811 0 -0.0114113 + 0.103507 0 -0.0114518 + 0.105321 0 -0.0114725 + 0.10713 0 -0.0114691 + 0.108808 0 -0.0114411 + 0.110579 0 -0.0113781 + 0.112302 0 -0.0112626 + 0.11382 0 -0.0111141 + 0.115346 0 -0.0112762 + 0.11714 0 -0.0114078 + 0.119066 0 -0.0114933 + 0.120972 0 -0.0115532 + 0.122742 0 -0.0115961 + 0.124611 0 -0.0116185 + 0.126457 0 -0.0116145 + 0.12817 0 -0.0115938 + 0.105384 0 -0.0132316 + 0.110716 0 -0.0131268 + 0.112608 0 -0.012929 + 0.113777 0 -0.0123264 + 0.114905 0 -0.0129363 + 0.116859 0 -0.0131428 + 0.119028 0 -0.0132271 + 0.124577 0 -0.0133094 + 0.103224 0 -0.00397825 + 0.105138 0 -0.00378865 + 0.106004 0 -0.00376124 + 0.101528 0 -0.00370832 + 0.102597 0 -0.00382751 + 0.101039 0 -0.00341169 + 0.10082 0 -0.0041226 + 0.103439 0 -0.00346991 + 0.10059 0 -0.0035272 + 0.101984 0 -0.0132038 + 0.103401 0 -0.0132218 + 0.107366 0 -0.0132251 + 0.108768 0 -0.0131983 + 0.121096 0 -0.0132702 + 0.122558 0 -0.0132939 + 0.12658 0 -0.0133066 + 0.128006 0 -0.0132964 + 0.104149 0 -0.00358053 + 0.102156 0 -0.00338396 + 0.102884 0 -0.00332436 + 0.101162 0 -0.00497468 + 0.102126 0 -0.00462898 + 0.101707 0 -0.00599016 + 0.085 0 0.0272727 + 0.085 0 0.0245455 + 0.085 0 0.0218182 + 0.085 0 0.0190909 + 0.085 0 0.0163636 + 0.085 0 0.0136364 + 0.085 0 0.0109091 + 0.085 0 0.00818182 + 0.085 0 0.00545455 + 0.085 0 0.00272727 + 0.0868219 0 0 + 0.0886438 0 0 + 0.0904656 0 0 + 0.0922304 0 0 + 0.0937169 0 0 + 0.0948407 0 0 + 0.0957179 0 0 + 0.0964255 0 0 + 0.0952885 0 0.03 + 0.0935737 0 0.03 + 0.091859 0 0.03 + 0.0901442 0 0.03 + 0.0884295 0 0.03 + 0.0867147 0 0.03 + 0.1 0 -0.00357452 + 0.1 0 -0.00428206 + 0.1 0 -0.00515926 + 0.1 0 -0.00628309 + 0.1 0 -0.00776959 + 0.1 0 -0.00953435 + 0.1 0 -0.0113562 + 0.1 0 -0.0131781 + 0.102727 0 -0.015 + 0.105455 0 -0.015 + 0.108182 0 -0.015 + 0.110909 0 -0.015 + 0.113636 0 -0.015 + 0.116364 0 -0.015 + 0.119091 0 -0.015 + 0.121818 0 -0.015 + 0.124545 0 -0.015 + 0.127273 0 -0.015 + 0.13 0 -0.0132853 + 0.13 0 -0.0115705 + 0.13 0 -0.00985577 + 0.13 0 -0.00814103 + 0.13 0 -0.00642629 + 0.13 0 -0.00471154 + 0.085 0 0.03 + 0.085 0 0 + 0.1 0 -0.015 + 0.13 0 -0.015 + + + 598 300 301 + 599 329 330 + 326 327 600 + 335 336 601 + 604 603 602 + 367 368 605 + 370 371 606 + 899 779 778 + 785 784 908 + 901 781 780 + 783 782 906 + 788 787 786 + 147 789 148 + 893 655 206 + 198 199 894 + 898 897 896 + 931 886 885 + 892 891 940 + 933 888 887 + 890 889 938 + 197 198 895 + + + 0 21 1 22 + 1 22 2 23 + 2 23 3 24 + 3 24 4 25 + 4 25 5 26 + 5 26 6 27 + 6 27 7 28 + 7 28 8 29 + 8 29 9 30 + 9 30 10 31 + 10 31 11 32 + 11 32 12 33 + 12 33 13 34 + 13 34 14 35 + 14 35 15 36 + 15 36 16 37 + 16 37 17 38 + 17 38 18 39 + 18 39 19 40 + 19 40 20 41 + 21 42 22 43 + 22 43 23 44 + 23 44 24 45 + 24 45 25 46 + 25 46 26 47 + 26 47 27 48 + 27 48 28 49 + 28 49 29 50 + 29 50 30 51 + 30 51 31 52 + 31 52 32 53 + 32 53 33 54 + 33 54 34 55 + 34 55 35 56 + 35 56 36 57 + 36 57 37 58 + 37 58 38 59 + 38 59 39 60 + 39 60 40 61 + 40 61 41 62 + 42 63 43 64 + 43 64 44 65 + 44 65 45 66 + 45 66 46 67 + 46 67 47 68 + 47 68 48 69 + 48 69 49 70 + 49 70 50 71 + 50 71 51 72 + 51 72 52 73 + 52 73 53 74 + 53 74 54 75 + 54 75 55 76 + 55 76 56 77 + 56 77 57 78 + 57 78 58 79 + 58 79 59 80 + 59 80 60 81 + 60 81 61 82 + 61 82 62 83 + 63 84 64 85 + 64 85 65 86 + 65 86 66 87 + 66 87 67 88 + 67 88 68 89 + 68 89 69 90 + 69 90 70 91 + 70 91 71 92 + 71 92 72 93 + 72 93 73 94 + 73 94 74 95 + 74 95 75 96 + 75 96 76 97 + 76 97 77 98 + 77 98 78 99 + 78 99 79 100 + 79 100 80 101 + 80 101 81 102 + 81 102 82 103 + 82 103 83 104 + 84 105 85 106 + 85 106 86 107 + 86 107 87 108 + 87 108 88 109 + 88 109 89 110 + 89 110 90 111 + 90 111 91 112 + 91 112 92 113 + 92 113 93 114 + 93 114 94 115 + 94 115 95 116 + 95 116 96 117 + 96 117 97 118 + 97 118 98 119 + 98 119 99 120 + 99 120 100 121 + 100 121 101 122 + 101 122 102 123 + 102 123 103 124 + 103 124 104 125 + 105 126 106 127 + 106 127 107 128 + 107 128 108 129 + 108 129 109 130 + 109 130 110 131 + 110 131 111 132 + 111 132 112 133 + 112 133 113 134 + 113 134 114 135 + 114 135 115 136 + 115 136 116 137 + 116 137 117 138 + 117 138 118 139 + 118 139 119 140 + 119 140 120 141 + 120 141 121 142 + 121 142 122 143 + 122 143 123 144 + 123 144 124 145 + 124 145 125 146 + 205 202 174 0 + 202 201 0 21 + 201 200 21 42 + 200 199 42 63 + 199 198 63 84 + 198 197 84 105 + 197 196 105 126 + 196 206 126 175 + 126 175 127 176 + 127 176 128 177 + 128 177 129 178 + 129 178 130 179 + 130 179 131 180 + 131 180 132 181 + 132 181 133 182 + 133 182 134 183 + 134 183 135 184 + 135 184 136 185 + 136 185 137 186 + 137 186 138 187 + 138 187 139 188 + 139 188 140 189 + 140 189 141 190 + 141 190 142 191 + 142 191 143 192 + 143 192 144 193 + 144 193 145 194 + 145 194 146 195 + 146 195 147 204 + 125 146 148 147 + 104 125 149 148 + 83 104 150 149 + 62 83 151 150 + 41 62 152 151 + 20 41 153 152 + 154 20 203 153 + 155 19 154 20 + 156 18 155 19 + 157 17 156 18 + 158 16 157 17 + 159 15 158 16 + 160 14 159 15 + 161 13 160 14 + 162 12 161 13 + 163 11 162 12 + 164 10 163 11 + 165 9 164 10 + 166 8 165 9 + 167 7 166 8 + 168 6 167 7 + 169 5 168 6 + 170 4 169 5 + 171 3 170 4 + 172 2 171 3 + 173 1 172 2 + 174 0 173 1 + 207 678 208 623 + 208 623 209 624 + 209 624 210 625 + 210 625 211 626 + 211 626 212 627 + 212 627 213 628 + 213 628 214 629 + 214 629 215 630 + 215 630 216 631 + 216 631 217 632 + 217 632 218 633 + 218 633 219 634 + 219 634 220 635 + 220 635 221 636 + 221 636 222 637 + 222 637 223 638 + 223 638 224 639 + 224 639 607 677 + 224 607 225 608 + 225 608 226 609 + 226 609 227 610 + 227 610 228 611 + 228 611 229 612 + 229 612 230 613 + 230 613 231 614 + 231 614 232 615 + 232 615 233 616 + 233 616 234 617 + 234 617 235 618 + 235 618 236 619 + 236 619 237 620 + 237 620 238 621 + 238 621 239 622 + 239 622 195 204 + 239 195 240 194 + 240 194 241 193 + 241 193 242 192 + 242 192 243 191 + 243 191 244 190 + 244 190 245 189 + 245 189 246 188 + 246 188 247 187 + 247 187 248 186 + 248 186 249 185 + 249 185 250 184 + 250 184 251 183 + 251 183 252 182 + 252 182 253 181 + 253 181 254 180 + 254 180 255 179 + 255 179 256 178 + 256 178 257 177 + 257 177 258 176 + 258 176 259 175 + 259 175 655 206 + 259 655 260 654 + 260 654 261 653 + 261 653 262 652 + 262 652 263 651 + 263 651 264 650 + 264 650 265 649 + 265 649 266 648 + 266 648 267 647 + 267 647 268 646 + 268 646 269 645 + 269 645 270 644 + 270 644 271 643 + 271 643 272 642 + 272 642 273 641 + 273 641 274 640 + 274 640 660 679 + 274 660 275 661 + 275 661 276 662 + 276 662 277 663 + 277 663 278 664 + 278 664 279 665 + 279 665 280 666 + 280 666 281 667 + 281 667 282 668 + 282 668 283 669 + 283 669 284 670 + 284 670 285 671 + 285 671 286 672 + 286 672 287 673 + 287 673 288 674 + 288 674 289 675 + 289 675 290 676 + 290 676 291 680 + 291 680 292 656 + 292 656 293 657 + 293 657 294 658 + 294 658 295 659 + 295 659 207 678 + 296 207 297 208 + 297 208 298 209 + 298 209 299 210 + 299 210 300 211 + 300 211 301 212 + 301 212 302 213 + 302 213 303 214 + 303 214 304 215 + 304 215 305 216 + 305 216 306 217 + 306 217 307 218 + 307 218 308 219 + 308 219 309 220 + 309 220 310 221 + 310 221 311 222 + 311 222 312 223 + 312 223 225 224 + 312 225 313 226 + 313 226 314 227 + 314 227 315 228 + 315 228 316 229 + 316 229 317 230 + 317 230 318 231 + 318 231 319 232 + 319 232 320 233 + 320 233 321 234 + 321 234 322 235 + 322 235 323 236 + 323 236 324 237 + 324 237 325 238 + 325 238 240 239 + 325 240 326 241 + 326 241 327 242 + 327 242 328 243 + 328 243 329 244 + 329 244 330 245 + 330 245 331 246 + 331 246 332 247 + 332 247 333 248 + 333 248 334 249 + 334 249 335 250 + 335 250 336 251 + 336 251 337 252 + 337 252 338 253 + 338 253 339 254 + 339 254 340 255 + 340 255 341 256 + 341 256 342 257 + 342 257 343 258 + 343 258 260 259 + 343 260 344 261 + 344 261 345 262 + 345 262 346 263 + 346 263 347 264 + 347 264 348 265 + 348 265 349 266 + 349 266 350 267 + 350 267 351 268 + 351 268 352 269 + 352 269 353 270 + 353 270 354 271 + 354 271 355 272 + 355 272 356 273 + 356 273 275 274 + 356 275 357 276 + 357 276 358 277 + 358 277 359 278 + 359 278 360 279 + 360 279 361 280 + 361 280 362 281 + 362 281 363 282 + 363 282 364 283 + 364 283 365 284 + 365 284 366 285 + 366 285 367 286 + 367 286 368 287 + 368 287 369 288 + 369 288 370 289 + 370 289 371 290 + 371 290 372 291 + 372 291 373 292 + 373 292 374 293 + 374 293 375 294 + 375 294 376 295 + 376 295 296 207 + 377 392 378 393 + 380 394 381 395 + 381 395 382 396 + 382 396 383 397 + 383 397 384 398 + 384 398 385 399 + 385 399 386 400 + 386 400 387 401 + 387 401 388 402 + 388 402 389 403 + 389 403 390 404 + 390 404 391 405 + 391 405 392 406 + 392 406 393 407 + 394 408 395 409 + 395 409 396 410 + 396 410 397 411 + 397 411 398 412 + 398 412 399 413 + 399 413 400 414 + 400 414 401 415 + 401 415 402 416 + 402 416 403 417 + 403 417 404 418 + 404 418 405 419 + 405 419 406 420 + 406 420 407 421 + 408 422 409 423 + 409 423 410 424 + 410 424 411 425 + 411 425 412 426 + 412 426 413 427 + 413 427 414 428 + 414 428 415 429 + 415 429 416 430 + 416 430 417 431 + 417 431 418 432 + 418 432 419 433 + 419 433 420 434 + 420 434 421 435 + 422 436 423 437 + 423 437 424 438 + 424 438 425 439 + 425 439 426 440 + 426 440 427 441 + 427 441 428 442 + 428 442 429 443 + 429 443 430 444 + 430 444 431 445 + 431 445 432 446 + 432 446 433 447 + 433 447 434 448 + 434 448 435 449 + 436 450 437 451 + 437 451 438 452 + 438 452 439 453 + 439 453 440 454 + 440 454 441 455 + 441 455 442 456 + 442 456 443 457 + 443 457 444 458 + 444 458 445 459 + 445 459 446 460 + 446 460 447 461 + 447 461 448 462 + 448 462 449 463 + 450 464 451 465 + 451 465 452 466 + 452 466 453 467 + 453 467 454 468 + 454 468 455 469 + 455 469 456 470 + 456 470 457 471 + 457 471 458 472 + 458 472 459 473 + 459 473 460 474 + 460 474 461 475 + 461 475 462 476 + 462 476 463 477 + 464 478 465 479 + 465 479 466 480 + 466 480 467 481 + 467 481 468 482 + 468 482 469 483 + 469 483 470 484 + 470 484 471 485 + 471 485 472 486 + 472 486 473 487 + 473 487 474 488 + 474 488 475 489 + 475 489 476 490 + 476 490 477 491 + 478 492 479 493 + 479 493 480 494 + 480 494 481 495 + 481 495 482 496 + 482 496 483 497 + 483 497 484 498 + 484 498 485 499 + 485 499 486 500 + 486 500 487 501 + 487 501 488 502 + 488 502 489 503 + 489 503 490 504 + 490 504 491 505 + 492 506 493 507 + 493 507 494 508 + 494 508 495 509 + 495 509 496 510 + 496 510 497 511 + 497 511 498 512 + 498 512 499 513 + 499 513 500 514 + 500 514 501 515 + 501 515 502 516 + 502 516 503 517 + 503 517 504 518 + 504 518 505 519 + 506 520 507 521 + 507 521 508 522 + 508 522 509 523 + 509 523 510 524 + 510 524 511 525 + 511 525 512 526 + 512 526 513 527 + 513 527 514 528 + 514 528 515 529 + 515 529 516 530 + 516 530 517 531 + 517 531 518 532 + 518 532 519 533 + 520 534 521 535 + 521 535 522 536 + 522 536 523 537 + 523 537 524 538 + 524 538 525 539 + 525 539 526 540 + 526 540 527 541 + 527 541 528 542 + 528 542 529 543 + 529 543 530 544 + 530 544 531 545 + 531 545 532 546 + 532 546 533 547 + 535 548 536 549 + 536 549 537 550 + 537 550 538 551 + 538 551 539 552 + 539 552 540 553 + 540 553 541 554 + 541 554 542 555 + 542 555 543 556 + 543 556 544 557 + 544 557 545 558 + 545 558 546 559 + 546 559 547 560 + 549 561 550 562 + 550 562 551 563 + 551 563 552 564 + 552 564 553 565 + 553 565 554 566 + 554 566 555 567 + 555 567 556 568 + 556 568 557 569 + 557 569 558 570 + 558 570 559 571 + 559 571 560 572 + 384 385 307 306 + 307 308 384 383 + 308 309 383 382 + 309 310 382 381 + 310 311 381 380 + 311 312 380 313 + 303 304 388 387 + 304 305 387 386 + 305 306 386 385 + 301 302 390 389 + 302 303 389 388 + 598 301 391 390 + 392 377 391 598 + 377 299 598 300 + 574 576 520 534 + 576 599 534 573 + 599 575 329 328 + 599 330 573 331 + 464 450 319 318 + 450 436 318 317 + 436 422 317 316 + 422 408 316 315 + 408 394 315 314 + 394 380 314 313 + 327 328 600 575 + 326 600 325 324 + 324 600 323 578 + 600 575 578 574 + 574 520 578 579 + 520 506 579 577 + 478 464 320 319 + 492 478 321 320 + 321 577 492 506 + 333 334 580 548 + 581 549 601 548 + 581 601 337 336 + 604 581 338 337 + 335 601 334 548 + 349 350 566 567 + 350 351 567 568 + 351 352 568 569 + 352 353 569 570 + 353 354 570 571 + 354 355 571 572 + 603 583 602 584 + 563 562 582 585 + 348 349 565 566 + 347 348 564 565 + 339 603 338 604 + 370 606 369 586 + 606 393 586 407 + 367 605 435 421 + 369 586 368 605 + 586 407 605 421 + 435 449 367 366 + 449 463 366 365 + 463 477 365 364 + 477 491 364 363 + 491 505 363 362 + 505 519 362 361 + 519 533 361 360 + 533 547 360 359 + 547 560 359 358 + 560 572 358 357 + 572 355 357 356 + 374 375 379 376 + 298 299 378 377 + 575 599 574 576 + 579 577 322 321 + 579 322 578 323 + 580 548 587 535 + 587 535 573 534 + 333 580 332 587 + 332 587 331 573 + 581 604 549 561 + 583 603 340 339 + 345 582 584 585 + 561 604 588 602 + 602 584 588 585 + 585 562 588 561 + 341 589 340 583 + 341 342 589 590 + 342 343 590 344 + 590 344 584 345 + 589 590 583 584 + 296 591 376 379 + 296 297 591 592 + 297 298 592 378 + 346 347 593 564 + 564 563 593 582 + 582 345 593 346 + 606 371 594 372 + 594 372 595 373 + 595 373 379 374 + 592 378 596 393 + 393 606 596 594 + 596 594 597 595 + 597 595 591 379 + 591 592 597 596 + 778 779 681 686 + 681 686 682 687 + 682 687 683 688 + 683 688 684 689 + 684 689 685 690 + 779 691 686 692 + 686 692 687 693 + 687 693 688 694 + 688 694 689 695 + 689 695 690 696 + 691 780 692 697 + 692 697 693 698 + 693 698 694 699 + 694 699 695 700 + 695 700 696 701 + 780 781 697 702 + 697 702 698 703 + 698 703 699 704 + 699 704 700 705 + 700 705 701 706 + 781 707 702 708 + 702 708 703 709 + 703 709 704 710 + 704 710 705 711 + 705 711 706 712 + 707 713 708 714 + 708 714 709 715 + 709 715 710 716 + 710 716 711 717 + 711 717 712 718 + 713 719 714 720 + 714 720 715 721 + 715 721 716 722 + 716 722 717 723 + 717 723 718 724 + 719 725 720 726 + 720 726 721 727 + 721 727 722 728 + 722 728 723 729 + 723 729 724 730 + 725 731 726 732 + 726 732 727 733 + 727 733 728 734 + 728 734 729 735 + 729 735 730 736 + 731 737 732 738 + 732 738 733 739 + 733 739 734 740 + 734 740 735 741 + 735 741 736 742 + 737 782 738 743 + 738 743 739 744 + 739 744 740 745 + 740 745 741 746 + 741 746 742 747 + 782 783 743 748 + 743 748 744 749 + 744 749 745 750 + 745 750 746 751 + 746 751 747 752 + 783 753 748 754 + 748 754 749 755 + 749 755 750 756 + 750 756 751 757 + 751 757 752 758 + 753 784 754 759 + 754 759 755 760 + 755 760 756 761 + 756 761 757 762 + 757 762 758 763 + 784 785 759 764 + 759 764 760 765 + 760 765 761 766 + 761 766 762 786 + 762 786 763 787 + 768 769 767 787 + 607 677 685 917 + 608 607 690 685 + 696 701 609 610 + 609 608 696 690 + 763 771 758 770 + 622 621 771 770 + 772 770 620 621 + 752 772 619 620 + 611 610 706 701 + 706 712 611 612 + 613 612 718 712 + 718 724 613 614 + 747 752 618 619 + 615 614 730 724 + 616 615 736 730 + 742 747 617 618 + 617 616 742 736 + 684 685 918 917 + 918 919 684 683 + 682 683 920 919 + 921 922 681 778 + 681 682 921 920 + 899 900 779 691 + 899 778 947 922 + 784 753 908 907 + 785 908 909 948 + 901 902 781 707 + 901 780 900 691 + 782 737 906 905 + 783 906 753 907 + 903 713 902 707 + 731 904 737 905 + 764 785 910 909 + 765 764 911 910 + 912 913 766 786 + 766 765 912 911 + 153 773 203 916 + 788 786 914 913 + 774 767 788 787 + 915 916 774 773 + 774 788 915 914 + 789 763 769 787 + 149 769 150 768 + 771 763 775 789 + 775 789 204 147 + 775 204 771 622 + 772 752 770 758 + 903 904 725 731 + 903 725 713 719 + 773 153 776 152 + 776 767 773 774 + 789 769 148 149 + 776 152 777 151 + 151 150 777 768 + 768 767 777 776 + 897 898 790 805 + 790 805 791 806 + 791 806 792 807 + 792 807 793 808 + 793 808 794 809 + 794 809 795 810 + 795 810 796 811 + 796 811 797 812 + 797 812 798 813 + 798 813 799 814 + 799 814 800 815 + 800 815 801 816 + 801 816 802 817 + 802 817 803 818 + 803 818 804 819 + 898 820 805 821 + 805 821 806 822 + 806 822 807 823 + 807 823 808 824 + 808 824 809 825 + 809 825 810 826 + 810 826 811 827 + 811 827 812 828 + 812 828 813 829 + 813 829 814 830 + 814 830 815 831 + 815 831 816 832 + 816 832 817 833 + 817 833 818 834 + 818 834 819 835 + 820 836 821 837 + 821 837 822 838 + 822 838 823 839 + 823 839 824 840 + 824 840 825 841 + 825 841 826 842 + 826 842 827 843 + 827 843 828 844 + 828 844 829 845 + 829 845 830 846 + 830 846 831 847 + 831 847 832 848 + 832 848 833 849 + 833 849 834 850 + 834 850 835 851 + 836 852 837 853 + 837 853 838 854 + 838 854 839 855 + 839 855 840 856 + 840 856 841 857 + 841 857 842 858 + 842 858 843 859 + 843 859 844 860 + 844 860 845 861 + 845 861 846 862 + 846 862 847 863 + 847 863 848 864 + 848 864 849 865 + 849 865 850 866 + 850 866 851 867 + 852 885 853 886 + 853 886 854 868 + 854 868 855 887 + 855 887 856 888 + 856 888 857 869 + 857 869 858 870 + 858 870 859 871 + 859 871 860 872 + 860 872 861 873 + 861 873 862 874 + 862 874 863 889 + 863 889 864 890 + 864 890 865 875 + 865 875 866 891 + 866 891 867 892 + 803 804 641 640 + 642 643 802 801 + 802 803 642 641 + 893 790 877 791 + 893 877 655 654 + 877 878 654 653 + 878 792 653 652 + 800 801 644 643 + 644 645 800 799 + 798 799 646 645 + 646 647 798 797 + 651 652 793 792 + 796 797 648 647 + 795 796 649 648 + 650 651 794 793 + 794 795 650 649 + 880 894 897 879 + 199 200 894 879 + 929 930 852 885 + 928 929 836 852 + 820 898 927 926 + 927 928 820 836 + 881 879 201 200 + 898 896 926 925 + 879 882 897 896 + 924 925 882 896 + 931 932 886 868 + 931 885 949 930 + 891 875 940 939 + 892 940 941 950 + 933 934 888 869 + 933 887 932 868 + 889 874 938 937 + 890 938 875 939 + 935 870 934 869 + 873 936 874 937 + 819 835 945 944 + 946 679 804 640 + 804 819 946 945 + 943 944 851 835 + 867 892 942 941 + 942 943 867 851 + 198 894 895 880 + 880 897 876 790 + 876 790 883 893 + 883 893 196 206 + 880 876 895 883 + 895 883 197 196 + 791 792 877 878 + 881 201 884 202 + 884 202 923 205 + 879 881 882 884 + 882 884 924 923 + 935 936 872 873 + 935 872 870 871 + + From 84f83b658efbc9fd57661987292df5e2c80280fa Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Thu, 6 Nov 2025 19:39:47 +0100 Subject: [PATCH 06/11] Fixed: IOP=2 option for SHLTRI and input to SHLTRI from Q4_INV. Tighten the comparison with zero tolerance (1e-8 is too loose). --- src/ASM/Ftn/dnorm3.f | 4 ++-- src/ASM/Ftn/q4_inv.f | 4 ++-- src/ASM/Ftn/shltri.f | 12 ++++++++---- 3 files changed, 12 insertions(+), 8 deletions(-) diff --git a/src/ASM/Ftn/dnorm3.f b/src/ASM/Ftn/dnorm3.f index 441969c69..55c5fa2da 100644 --- a/src/ASM/Ftn/dnorm3.f +++ b/src/ASM/Ftn/dnorm3.f @@ -3,8 +3,8 @@ DOUBLE PRECISION FUNCTION DNORM3 (V) C C DNORM3 : Calculate |v| and 1/|v| * {v} C - DOUBLE PRECISION ZERO , EPS , ONE - PARAMETER ( ZERO = 0.0D0, EPS = 1.0D-8, ONE = 1.0D0 ) + DOUBLE PRECISION ZERO , EPS , ONE + PARAMETER ( ZERO = 0.0D0, EPS = 1.0D-15, ONE = 1.0D0 ) C DV = SQRT(V(1)*V(1) + V(2)*V(2) + V(3)*V(3)) IF (DV .LT. EPS) THEN diff --git a/src/ASM/Ftn/q4_inv.f b/src/ASM/Ftn/q4_inv.f index 77dcb674c..ee46ebb02 100644 --- a/src/ASM/Ftn/q4_inv.f +++ b/src/ASM/Ftn/q4_inv.f @@ -123,8 +123,8 @@ SUBROUTINE Q4_INV (NDIM,EPSO,EPSZ,X,TENC,XI,IPSW,IWR,IERR) C to the global XY-plane, need C to calculate local coordinates DO 2000 j = 1,NDIM - X1(j) = TENC(j,3) - TENC(j,1) - X2(j) = TENC(j,4) - TENC(j,2) + X1(j) = TENC(j,2) - TENC(j,1) + X2(j) = TENC(j,3) - TENC(j,1) 2000 CONTINUE CALL SHLTRI (2,X1,X2,TGL,IPSW,IWR,IERRL) IF (IERRL .LT. 0) GOTO 7100 diff --git a/src/ASM/Ftn/shltri.f b/src/ASM/Ftn/shltri.f index 6b5c8b38c..ceca4dc0f 100644 --- a/src/ASM/Ftn/shltri.f +++ b/src/ASM/Ftn/shltri.f @@ -81,7 +81,7 @@ SUBROUTINE SHLTRI (IOP,VX,V3,TGL,IPSW,IWR,IERR) IF (IPSW .GE. 6) WRITE(IWR,9900) IF (IPSW .GE. 7) THEN WRITE(IWR,9920) IOP,V3 - IF (IOP .EQ. 1) WRITE(IWR,9921) VX + IF (IOP .EQ. 1 .OR. IOP .EQ. 2) WRITE(IWR,9921) VX ENDIF C C @@ -106,13 +106,17 @@ SUBROUTINE SHLTRI (IOP,VX,V3,TGL,IPSW,IWR,IERR) C ELSE IF (IOP .EQ. 2) THEN C -C * Define V2 as V1 x V3, -C then V3 = V1 x V2 + V1(1) = VX(1) + V1(2) = VX(2) + V1(3) = VX(3) +C +C * Define V3 as V1 x V3, +C then V2 = V3 x V1 CALL DCROS3 (V1,V3,V2) V3(1) = V2(1) V3(2) = V2(2) V3(3) = V2(3) - CALL DCROS3 (V1,V2,V3) + CALL DCROS3 (V3,V1,V2) IF (DNORM3(V3) .LE. ZERO) GOTO 7200 C ELSE IF (ABS(V3(2)) .GT. EPS .OR. ABS(V3(3)) .GT. EPS) THEN From 5f138cca321bfbb9a1524c21dbfb3debbe90f9ef Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Sat, 8 Nov 2025 11:01:04 +0100 Subject: [PATCH 07/11] 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 2096b8df93d91d79b207efedc1b4052624a31c6d Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Sun, 9 Nov 2025 07:19:37 +0100 Subject: [PATCH 08/11] =?UTF-8?q?Added:=20Another=20(non-virtual)=C2=A0eva?= =?UTF-8?q?lSolution()=20method=20evaluating=20a=20nodal=20field=20at=20sp?= =?UTF-8?q?ecified=20points=20defined=20by=20element=20index=20and=20local?= =?UTF-8?q?=20xi,eta=20parameters=20within=20those=20element.=20Added:=20p?= =?UTF-8?q?rivate=20method=20evalPt()=20evaluating=20a=20nodal=20field=20a?= =?UTF-8?q?t=20a=20point.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/ASM/ASMs2DLag.C | 73 ++++++++++++++++++++++++++++++++++----------- src/ASM/ASMs2DLag.h | 18 +++++++++++ src/ASM/ASMs3DLag.C | 47 +++++++++++++++++++---------- src/ASM/ASMs3DLag.h | 4 +++ 4 files changed, 109 insertions(+), 33 deletions(-) diff --git a/src/ASM/ASMs2DLag.C b/src/ASM/ASMs2DLag.C index 0fb25a6d4..adbaffa91 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,11 +742,9 @@ 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)) @@ -757,20 +756,12 @@ bool ASMs2DLag::evalSolution (Matrix& sField, const Vector& locSol, for (size_t i = 0; i < ni; i++, ip++) { 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++; - - for (size_t a = 1; a <= nen; a++) - elmSol.fillColumn(a, locSol.ptr() + nCmp*MNPC[iel-1][a-1]); + else // evaluate at element centers + --iel; + if (!this->evalSolPt(iel,xi,eta,nCmp,locSol,val,N)) + return false; - elmSol.multiply(N,val); sField.fillColumn(ip,val); } @@ -778,6 +769,54 @@ bool ASMs2DLag::evalSolution (Matrix& sField, const Vector& locSol, } +bool ASMs2DLag::evalSolution (Matrix& sField, const Vector& locSol, + const IntVec& elms, const RealArray* lpar) const +{ + const size_t nCmp = locSol.size() / this->getNoProjectionNodes(); + const size_t nPts = lpar ? lpar[0].size() : elms.size(); + const size_t iPts = sField.cols(); + sField.resize(nCmp,iPts+nPts); + + RealArray N, ptSol; + for (size_t i = 0; i < nPts; i++) + { + double xi = lpar ? lpar[0][i] : 0.0; + double eta = lpar ? lpar[1][i] : 0.0; + if (!this->evalSolPt(elms[i],xi,eta,nCmp,locSol,ptSol,N)) + return false; + + sField.fillColumn(iPts+i+1,ptSol); + } + + return true; +} + + +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 (MNPC[--iel].size() == 3) + N = { xi, eta, 1.0-xi-eta }; // special for linear triangles + else if (!Lagrange::computeBasis(N,p1,xi,p2,eta)) + return false; + + 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; +} + + bool ASMs2DLag::evalSolution (Matrix& sField, const IntegrandBase& integrand, const int*, char) const { @@ -997,7 +1036,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..b049d00e8 100644 --- a/src/ASM/ASMs2DLag.h +++ b/src/ASM/ASMs2DLag.h @@ -233,6 +233,20 @@ class ASMs2DLag : public ASMs2D, protected ASMLagBase const RealArray* gpar, bool regular, int, int) const; + //! \brief Evaluates the primary solution field at the given points. + //! \param sField Solution field + //! \param[in] locSol Solution vector local to current patch + //! \param[in] elms List of elements to evaluate at + //! \param[in] lpar Local parameter values of the result sampling points + //! + //! \details If \a lpar is null, the solution is evaluated at the center of + //! each element. Otherwise, it is assumed to contain the local ξ and η + //! parameters of the evaluation points, w.r.t. the provided elements. + //! \note The resuls are appended as extra columns to output matrix \a sField. + //! Therefore, any existing results is preserved. + bool evalSolution(Matrix& sField, const Vector& locSol, + const IntVec& elms, const RealArray* lpar = nullptr) const; + //! \brief Evaluates the secondary solution field at all visualization points. //! \param[out] sField Solution field //! \param[in] integrand Object with problem-specific data and methods @@ -289,6 +303,10 @@ class ASMs2DLag : public ASMs2D, protected ASMLagBase //! \brief Returns matrix of nodal point correspondance for a structured grid. static void createMNPC(size_t nx, size_t ny, int p1, int p2, IntMat& MNPC); + //! \brief Evaluates a nodal solution field at specified point in an element. + bool evalSolPt(int iel, double xi, double eta, size_t nCmp, + const Vector& pchSol, RealArray& ptSol, RealArray& N) const; + protected: size_t nx; //!< Number of nodes in first parameter direction size_t ny; //!< Number of nodes in second parameter direction diff --git a/src/ASM/ASMs3DLag.C b/src/ASM/ASMs3DLag.C index c210bf8c0..4809e0c80 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,11 +1001,9 @@ 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)) @@ -1018,21 +1016,13 @@ bool ASMs3DLag::evalSolution (Matrix& sField, const Vector& locSol, for (size_t i = 0; i < ni; i++, ip++) { 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++; - - for (size_t a = 1; a <= nen; a++) - elmSol.fillColumn(a, locSol.ptr() + nCmp*MNPC[iel-1][a-1]); + else // evaluate at element centers + --iel; + if (!this->evalSolPt(iel,xi,eta,zeta,nCmp,locSol,val,N)) + return false; - elmSol.multiply(N,val); sField.fillColumn(ip,val); } @@ -1040,6 +1030,31 @@ bool ASMs3DLag::evalSolution (Matrix& sField, const Vector& locSol, } +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; +} + + bool ASMs3DLag::evalSolution (Matrix& sField, const IntegrandBase& integrand, const int*, char) const { diff --git a/src/ASM/ASMs3DLag.h b/src/ASM/ASMs3DLag.h index c49e5695e..f0ae72348 100644 --- a/src/ASM/ASMs3DLag.h +++ b/src/ASM/ASMs3DLag.h @@ -312,6 +312,10 @@ class ASMs3DLag : public ASMs3D, protected ASMLagBase static void createMNPC(size_t nx, size_t ny, size_t nz, int p1, int p2, int p3, IntMat& MNPC); + //! \brief Evaluates a nodal solution field at specified point in an element. + bool evalSolPt(int iel, double xi, double eta, double zeta, size_t nCmp, + const Vector& pchSol, RealArray& ptSol, RealArray& N) const; + protected: size_t nx; //!< Number of nodes in first parameter direction size_t ny; //!< Number of nodes in second parameter direction From 449ff2a2593aec5f54b16a45edb3dbb800e930a8 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Sun, 9 Nov 2025 07:25:16 +0100 Subject: [PATCH 09/11] Added: RealArray* lpar as argument to new evalSolution() method. Added: virtual method ASMbase::printElmInfo(). --- src/ASM/ASMbase.C | 2 +- src/ASM/ASMbase.h | 8 +++++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/ASM/ASMbase.C b/src/ASM/ASMbase.C index cd59c4563..225daf07f 100644 --- a/src/ASM/ASMbase.C +++ b/src/ASM/ASMbase.C @@ -1785,7 +1785,7 @@ bool ASMbase::evalSolution (Matrix&, const IntegrandBase&, bool ASMbase::evalSolution (Matrix& sField, const IntegrandBase& integrand, - const IntVec& elements) const + const IntVec& elements, const RealArray*) const { if (elements.empty()) return this->evalSolution(sField,integrand,nullptr,true); diff --git a/src/ASM/ASMbase.h b/src/ASM/ASMbase.h index 9e0512aa7..ab3687b97 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); @@ -731,13 +733,17 @@ class ASMbase //! \param[out] sField Solution field //! \param[in] integrand Object with problem-specific data and methods //! \param[in] elements List of elements to evaluate for (all if empty) + //! \param[in] lpar Local parameter values of the result sampling points //! //! \details The secondary solution is derived from the primary solution, //! which is assumed to be stored within the \a integrand for current patch. //! This method evaluates the secondary solution at the center of the //! provided list of \a elements, or for all if the list is empty. + //! If \a lpar is specified, it is assumed to contain the local parameters + //! (w.r.t. the elements) of the evaluation points virtual bool evalSolution(Matrix& sField, const IntegrandBase& integrand, - const IntVec& elements) const; + const IntVec& elements, + const RealArray* lpar = nullptr) const; //! \brief Projects the secondary solution using a (discrete) global L2-fit. //! \param[out] sField Secondary solution field control point values From df6fc735d08636affab627f3381d8258e05200f3 Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Sun, 9 Nov 2025 07:40:59 +0100 Subject: [PATCH 10/11] Fixed: Missing argument to NextSiblingElement() calls. Added: Evaluation of 2ndary solution at specified Cartesian points. --- src/SIM/SIM1D.C | 2 +- src/SIM/SIM2D.C | 2 +- src/SIM/SIM3D.C | 2 +- src/SIM/SIMinput.C | 2 +- src/SIM/SIMoutput.C | 137 ++++++++++++++++++++++++++++++-------------- src/SIM/TimeStep.C | 2 +- 6 files changed, 100 insertions(+), 47 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 a5988d162..e9eb1ae57 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); @@ -138,10 +137,13 @@ bool SIMoutput::parseOutputTag (const tinyxml2::XMLElement* elem) points.back().second.emplace_back(newPoint); }; + Vec3 X0; // Optional coordinate shift for Cartesian coordinates + utl::getAttribute(elem,"shift",X0); + // 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) @@ -150,21 +152,31 @@ bool SIMoutput::parseOutputTag (const tinyxml2::XMLElement* elem) IFEM::cout <<"\tPoint "<< i <<": P"<< thePoint.patch; if (utl::getAttribute(point,"node",thePoint.inod) && thePoint.inod > 0) IFEM::cout <<" node = "<< thePoint.inod; + else if (point->FirstChild() && + utl::parseVec(thePoint.X,point->FirstChild()->Value())) + { + // Cartesian coordinates of a spatial point is specified, + thePoint.iel = -1; // need to search for the matching element/node + thePoint.X += X0; + IFEM::cout <<" X = "<< thePoint.X; + } 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 +218,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 +244,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 +421,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) @@ -426,23 +434,52 @@ void SIMoutput::preprocessResPtGroup (std::string& ptFile, ResPointVec& points) nX = 0; } + // Check if we have general spatial points input (Cartesian coordinates) + using PatchPoints = std::pair,Vec3Vec>; + std::map pchPoints; + if (pointMap.empty()) + for (ResultPoint& pt : points) + if (pt.iel < 0 || pt.inod < 0) + if (ASMbase* pch = this->getPatch(pt.patch,true); pch && !pch->empty()) + { + pchPoints[pch].second.emplace_back(pt.X); + pt.iel = -pchPoints[pch].second.size(); + } + + // Find element and/or local parameters of the spatial points + for (std::pair& pch : pchPoints) + if (!pch.first->findPoints(pch.second.second,pch.second.first)) + IFEM::cout <<" ** Failed to locate some result points."<< std::endl; + int iPoint = 0; size_t iX = 0, iY = 0; 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); - else if (pit->inod < 0) // A spatial point is specified - pointIsOK = fabs(pch->findPoint(pit->X,pit->u)) < 1.0e-3; - + else if (pit->iel < 0 && !pchPoints.empty()) + { + // A spatial point is specified, use the point searching results + size_t iP = -pit->iel-1; + std::map::const_iterator it = pchPoints.find(pch); + if (it != pchPoints.end() && iP < it->second.first.size()) + if ((pointIsOK = fabs(it->second.first[iP].dist) < 1.0e-3)) + { + pit->npar = pch->getNoParamDim(); + pit->X = it->second.second[iP]; + pit->iel = it->second.first[iP].iel; + memcpy(pit->u,it->second.first[iP].u,sizeof(double)*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 +492,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 +533,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; @@ -509,6 +558,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); } } @@ -539,7 +590,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; } @@ -561,8 +612,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; @@ -1796,11 +1846,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()) @@ -2021,7 +2071,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()) @@ -2124,6 +2174,8 @@ bool SIMoutput::evalResults (const Vectors& psol, const ResPointVec& gPoints, else if (pt.iel > 0) { elms.push_back(pt.iel); + for (short int k = 0; k < pt.npar; k++) + params[k].push_back(pt.u[k]); } else if (pt.npar < 0) { @@ -2226,7 +2278,8 @@ bool SIMoutput::evalResults (const Vectors& psol, const ResPointVec& gPoints, if (elms.empty()) ok = patch->evalSolution(tmp,*myProblem,params.data(),points.empty()); else // evaluate at specified elements - ok = patch->evalSolution(tmp,*myProblem,elms); + ok = patch->evalSolution(tmp,*myProblem,elms, + params.front().empty() ? nullptr : params.data()); if ((ok &= augment(sol2,tmp)) && getNames) for (size_t i = 0; i < myProblem->getNoFields(2); i++) 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 98382aa987d7e2c966fa12dfead01982fa2a653b Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Fri, 14 Nov 2025 14:26:16 +0100 Subject: [PATCH 11/11] Fixed: Skip 0D/1D elements when evaluating solution at element centers. Changed: Override evalSolPt() in ASMu2DLag with the triangle support. --- src/ASM/ASMs2DLag.C | 14 ++++++++------ src/ASM/ASMs2DLag.h | 7 ++++--- src/ASM/ASMu2DLag.C | 20 ++++++++++++++++++++ src/ASM/ASMu2DLag.h | 5 +++++ 4 files changed, 37 insertions(+), 9 deletions(-) diff --git a/src/ASM/ASMs2DLag.C b/src/ASM/ASMs2DLag.C index adbaffa91..4d7a3f045 100644 --- a/src/ASM/ASMs2DLag.C +++ b/src/ASM/ASMs2DLag.C @@ -750,21 +750,23 @@ bool ASMs2DLag::evalSolution (Matrix& sField, const Vector& locSol, 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); else // evaluate at element centers --iel; + if (!this->evalSolPt(iel,xi,eta,nCmp,locSol,val,N)) return false; - - sField.fillColumn(ip,val); + else if (!val.empty()) // skip elements with no results + sField.fillColumn(++ip,val); } + sField.resize(nCmp,ip); return true; } @@ -800,10 +802,10 @@ bool ASMs2DLag::evalSolPt (int iel, double xi, double eta, size_t nCmp, iel = -iel-1; else if (iel < 1 || iel > static_cast(nel)) return false; - else if (MNPC[--iel].size() == 3) - N = { xi, eta, 1.0-xi-eta }; // special for linear triangles 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++) diff --git a/src/ASM/ASMs2DLag.h b/src/ASM/ASMs2DLag.h index b049d00e8..76b9c7d70 100644 --- a/src/ASM/ASMs2DLag.h +++ b/src/ASM/ASMs2DLag.h @@ -303,11 +303,12 @@ class ASMs2DLag : public ASMs2D, protected ASMLagBase //! \brief Returns matrix of nodal point correspondance for a structured grid. 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. - bool evalSolPt(int iel, double xi, double eta, size_t nCmp, - const Vector& pchSol, RealArray& ptSol, RealArray& N) const; + virtual bool evalSolPt(int iel, double xi, double eta, size_t nCmp, + const Vector& pchSol, RealArray& ptSol, + RealArray& N) const; -protected: 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/ASMu2DLag.C b/src/ASM/ASMu2DLag.C index ec487651b..eb5e079ad 100644 --- a/src/ASM/ASMu2DLag.C +++ b/src/ASM/ASMu2DLag.C @@ -608,6 +608,26 @@ bool ASMu2DLag::integrate (Integrand& integrand, } +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 3fb75abb2..69cca9f4e 100644 --- a/src/ASM/ASMu2DLag.h +++ b/src/ASM/ASMu2DLag.h @@ -147,6 +147,11 @@ class ASMu2DLag : public ASMs2DLag bool writeXML(const char* fname) const; protected: + //! \brief Evaluates a nodal solution field at specified point in an element. + 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);