Skip to content

Commit

Permalink
Fix indexing bug. (#139)
Browse files Browse the repository at this point in the history 
* 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 <[email protected]>
  • Loading branch information
@ktbolt @MatteoSalvador
ktbolt and MatteoSalvador authored Nov 16, 2023
1 parent dcc7df9 commit c26b7c6
Show file tree
Hide file tree
Showing 5 changed files with 40 additions and 29 deletions.
30 changes: 20 additions & 10 deletions Code/Source/svFSI/all_fun.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<double>& s, int l, int u, bool pFlag)
Expand All @@ -305,6 +307,7 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, int dId, const Array<do
#ifdef debug_integ_v
DebugMsg dmsg(__func__, com_mod.cm.idcm());
dmsg.banner();
dmsg << "vInteg " << " ";
dmsg << "dId: " << dId;
dmsg << "l: " << l;
dmsg << "pFlag: " << pFlag;
Expand Down Expand Up @@ -544,6 +547,7 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co
#ifdef debug_integ_s
DebugMsg dmsg(__func__, com_mod.cm.idcm());
dmsg.banner();
dmsg << "IntegS " << " ";
dmsg << "lFa.iM: " << lFa.iM+1;
dmsg << "lFa.name: " << lFa.name;
dmsg << "lFa.eType: " << lFa.eType;
Expand Down Expand Up @@ -687,6 +691,8 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co
#ifdef debug_integ_V
DebugMsg dmsg(__func__, com_mod.cm.idcm());
dmsg.banner();
dmsg << "IntegV " << " ";
dmsg << "lFa.name: " << lFa.name;
#endif

auto& cm = com_mod.cm;
Expand Down Expand Up @@ -775,18 +781,20 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co

/// @brief This routine integrate s(l:u,:) over the surface faId.
///
/// Note that 'l' seems to be a length and 'uo' an offset. 'l' should never be 0.
/// Note that 'l' and 'u' should be 0-based and are used to index into 's'.
///
/// Reproduces 'FUNCTION IntegG(lFa, s, l, uo, THflag)'.
//
double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Array<double>& s, const int l, int uo, bool THflag)
double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa,
const Array<double>& s, const int l, std::optional<int> uo, bool THflag)
{
using namespace consts;

#define n_debug_integ_g
#ifdef debug_integ_g
DebugMsg dmsg(__func__, com_mod.cm.idcm());
dmsg.banner();
dmsg << "IntegG " << " ";
dmsg << "l: " << l;
dmsg << "uo: " << uo;
#endif
Expand All @@ -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;
Expand All @@ -819,16 +831,16 @@ double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, co
if (u-l+1 == nsd) {
Array<double> 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);

} else if (l == u) {
Vector<double> 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 {
Expand Down Expand Up @@ -1117,8 +1129,6 @@ void set_dmn_id(mshType& mesh, const int iDmn, const int ifirst, const int ilast
mesh.eId = Vector<int>(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;
Expand Down
3 changes: 2 additions & 1 deletion Code/Source/svFSI/all_fun.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@

#include "consts.h"

#include <optional>
#include <string>

namespace all_fun {
Expand All @@ -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<double>& s,
const int l, int uo=-1, bool THflag=false);
const int l, std::optional<int> uo=std::nullopt, bool THflag=false);

double integ(const ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const Array<double>& s);

Expand Down
2 changes: 1 addition & 1 deletion Code/Source/svFSI/baf_ini.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<double> 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)
Expand Down
22 changes: 11 additions & 11 deletions Code/Source/svFSI/set_bc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
}
Expand All @@ -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.
Expand Down
12 changes: 6 additions & 6 deletions Code/Source/svFSI/txt.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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");
}
Expand All @@ -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);
Expand Down

0 comments on commit c26b7c6

Please sign in to comment.