From 422c93696576d36dd664f49056508307c74982e6 Mon Sep 17 00:00:00 2001 From: Dave Parker Date: Wed, 15 Nov 2023 18:32:16 -0800 Subject: [PATCH] Fix indexing bug. (#139) * Fix indexing bug. * Fix indexing problem (again). * Modify the calls to all_fun::integ() functions to pass in an array's 0-based index range. --------- Co-authored-by: Matteo Salvador <44835004+MatteoSalvador@users.noreply.github.com> --- Code/Source/svFSI/all_fun.cpp | 30 ++++++++++++++++++++---------- Code/Source/svFSI/all_fun.h | 3 ++- Code/Source/svFSI/baf_ini.cpp | 2 +- Code/Source/svFSI/set_bc.cpp | 22 +++++++++++----------- Code/Source/svFSI/txt.cpp | 12 ++++++------ 5 files changed, 40 insertions(+), 29 deletions(-) diff --git a/Code/Source/svFSI/all_fun.cpp b/Code/Source/svFSI/all_fun.cpp index cd62e39..acb3b96 100644 --- a/Code/Source/svFSI/all_fun.cpp +++ b/Code/Source/svFSI/all_fun.cpp @@ -295,6 +295,8 @@ global(const ComMod& com_mod, const CmMod& cm_mod, const mshType& lM, const Arra /// @brief This routine integrate an equation over a particular domain /// +/// Note that 'l' and 'u' should be 0-based and are used to index into 's'. +// /// Replicates 'FUNCTION vInteg(dId, s, l, u, pFlag)' defined in ALLFUN.f. // double integ(const ComMod& com_mod, const CmMod& cm_mod, int dId, const Array& s, int l, int u, bool pFlag) @@ -305,6 +307,7 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, int dId, const Array& s, const int l, int uo, bool THflag) +double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, + const Array& s, const int l, std::optional uo, bool THflag) { using namespace consts; @@ -787,6 +794,7 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co #ifdef debug_integ_g DebugMsg dmsg(__func__, com_mod.cm.idcm()); dmsg.banner(); + dmsg << "IntegG " << " "; dmsg << "l: " << l; dmsg << "uo: " << uo; #endif @@ -796,9 +804,13 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co int insd = nsd - 1; int tnNo = com_mod.tnNo; - int u = l; - if (uo != -1) { - u = uo; + // Set u if uo is given. + int u{0}; + + if (uo.has_value()) { + u = uo.value(); + } else { + u = l; } bool flag = THflag; @@ -819,8 +831,8 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co if (u-l+1 == nsd) { Array vec(nsd,nNo); for (int a = 0; a < nNo; a++) { - for (int i = 0; i < nsd; i++) { - vec(i,a) = s(i+l-1,a); + for (int i = l, n = 0; i <= u; i++, n++) { + vec(n,a) = s(i,a); } } result = integ(com_mod, cm_mod, lFa, vec); @@ -828,7 +840,7 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co } else if (l == u) { Vector sclr(nNo); for (int a = 0; a < nNo; a++) { - sclr(a) = s(l-1,a); + sclr(a) = s(l,a); } result = integ(com_mod, cm_mod, lFa, sclr, flag); } else { @@ -1117,8 +1129,6 @@ void set_dmn_id(mshType& mesh, const int iDmn, const int ifirst, const int ilast mesh.eId = Vector(mesh.gnEl); } - //std::cout << "[set_dmn_id] first: " << first << std::endl; - // Set the iDimn'th bit for each element ID. for (int e = first; e <= last; e++) { mesh.eId[e] |= 1UL << iDmn; diff --git a/Code/Source/svFSI/all_fun.h b/Code/Source/svFSI/all_fun.h index 3eba86e..9830705 100644 --- a/Code/Source/svFSI/all_fun.h +++ b/Code/Source/svFSI/all_fun.h @@ -37,6 +37,7 @@ #include "consts.h" +#include #include namespace all_fun { @@ -61,7 +62,7 @@ namespace all_fun { bool pFlag=false); double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Array& s, - const int l, int uo=-1, bool THflag=false); + const int l, std::optional uo=std::nullopt, bool THflag=false); double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Array& s); diff --git a/Code/Source/svFSI/baf_ini.cpp b/Code/Source/svFSI/baf_ini.cpp index a960d66..3c843b2 100644 --- a/Code/Source/svFSI/baf_ini.cpp +++ b/Code/Source/svFSI/baf_ini.cpp @@ -305,7 +305,7 @@ void bc_ini(const ComMod& com_mod, const CmMod& cm_mod, bcType& lBc, faceType& l } else if (btest(lBc.bType, iBC_para)) { Vector center(3); for (int i = 0; i < nsd; i++) { - center(i) = all_fun::integ(com_mod, cm_mod, lFa, com_mod.x, i+1) / lFa.area; + center(i) = all_fun::integ(com_mod, cm_mod, lFa, com_mod.x, i) / lFa.area; } // gNodes is one if a node located on the boundary (beside iFa) diff --git a/Code/Source/svFSI/set_bc.cpp b/Code/Source/svFSI/set_bc.cpp index 8b793ca..1699454 100644 --- a/Code/Source/svFSI/set_bc.cpp +++ b/Code/Source/svFSI/set_bc.cpp @@ -103,8 +103,8 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod) auto& fa = com_mod.msh[iM].fa[iFa]; if (utils::btest(bc.bType, iBC_Neu)) { - cplBC.fa[ptr].Qo = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yo, 1, nsd); - cplBC.fa[ptr].Qn = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yn, 1, nsd); + cplBC.fa[ptr].Qo = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yo, 0, nsd-1); + cplBC.fa[ptr].Qn = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yn, 0, nsd-1); cplBC.fa[ptr].Po = 0.0; cplBC.fa[ptr].Pn = 0.0; #ifdef debug_calc_der_cpl_bc @@ -115,8 +115,8 @@ void calc_der_cpl_bc(ComMod& com_mod, const CmMod& cm_mod) } else if (utils::btest(bc.bType, iBC_Dir)) { double area = fa.area; - cplBC.fa[ptr].Po = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yo, nsd+1) / area; - cplBC.fa[ptr].Pn = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yn, nsd+1) / area; + cplBC.fa[ptr].Po = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yo, nsd) / area; + cplBC.fa[ptr].Pn = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yn, nsd) / area; cplBC.fa[ptr].Qo = 0.0; cplBC.fa[ptr].Qn = 0.0; #ifdef debug_calc_der_cpl_bc @@ -465,8 +465,8 @@ void rcr_init(ComMod& com_mod, const CmMod& cm_mod) if (cplBC.initRCR) { auto& fa = com_mod.msh[iM].fa[iFa]; double area = fa.area; - double Qo = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yo, 1, nsd); - double Po = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yo, nsd+1) / area; + double Qo = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yo, 0, nsd-1); + double Po = all_fun::integ(com_mod, cm_mod, fa, com_mod.Yo, nsd) / area; cplBC.xo[ptr] = Po - (Qo * cplBC.fa[ptr].RCR.Rp); } else { cplBC.xo[ptr] = cplBC.fa[ptr].RCR.Xo; @@ -592,16 +592,17 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod) } } + if (ptr != -1) { if (utils::btest(bc.bType,iBC_Neu)) { - cplBC.fa[ptr].Qo = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yo, 1, nsd); - cplBC.fa[ptr].Qn = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, 1, nsd); + cplBC.fa[ptr].Qo = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yo, 0, nsd-1); + cplBC.fa[ptr].Qn = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, 0, nsd-1); cplBC.fa[ptr].Po = 0.0; cplBC.fa[ptr].Pn = 0.0; } else if (utils::btest(bc.bType,iBC_Dir)) { double area = com_mod.msh[iM].fa[iFa].area; - cplBC.fa[ptr].Po = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yo, nsd+1) / area; - cplBC.fa[ptr].Pn = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, nsd+1) / area; + cplBC.fa[ptr].Po = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yo, nsd) / area; + cplBC.fa[ptr].Pn = all_fun::integ(com_mod, cm_mod, com_mod.msh[iM].fa[iFa], Yn, nsd) / area; cplBC.fa[ptr].Qo = 0.0; cplBC.fa[ptr].Qn = 0.0; } @@ -623,7 +624,6 @@ void set_bc_cpl(ComMod& com_mod, CmMod& cm_mod) bc.g = cplBC.fa[ptr].y; } } - } /// @brief Apply Dirichlet BCs strongly. diff --git a/Code/Source/svFSI/txt.cpp b/Code/Source/svFSI/txt.cpp index b2f48c7..1bfd972 100644 --- a/Code/Source/svFSI/txt.cpp +++ b/Code/Source/svFSI/txt.cpp @@ -477,17 +477,17 @@ void wtxt(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const int m, if (m == 1) { if (div) { tmp = fa.area; - tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 1) / tmp; + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0) / tmp; } else { if (pFlag && lTH) { - tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 1, true); + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, true); } else { - tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 1); + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0); } } } else if (m == nsd) { - tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 1, m); + tmp = all_fun::integ(com_mod, cm_mod, fa, tmpV, 0, m-1); } else { throw std::runtime_error("WTXT only accepts 1 and nsd"); } @@ -505,9 +505,9 @@ void wtxt(const ComMod& com_mod, CmMod& cm_mod, const eqType& lEq, const int m, if (div) { tmp = dmn.v; - tmp = all_fun::integ(com_mod, cm_mod, dmn.Id, tmpV, 1, m) / tmp; + tmp = all_fun::integ(com_mod, cm_mod, dmn.Id, tmpV, 0, m-1) / tmp; } else { - tmp = all_fun::integ(com_mod, cm_mod, dmn.Id, tmpV, 1, m, pFlag); + tmp = all_fun::integ(com_mod, cm_mod, dmn.Id, tmpV, 0, m-1, pFlag); } if (com_mod.cm.mas(cm_mod)) { fprintf(fp, " %.10e ", tmp);