Skip to content

Commit

Permalink
Split up construct_dsolid
Browse files Browse the repository at this point in the history
  • Loading branch information
@mrp089
mrp089 committed Dec 14, 2023
1 parent 5e604f8 commit a445e1f
Show file tree
Hide file tree
Showing 2 changed files with 135 additions and 113 deletions.
243 changes: 131 additions & 112 deletions Code/Source/svFSI/sv_struct.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -200,9 +200,67 @@ void b_struct_3d(const ComMod& com_mod, const int eNoN, const double w, const Ve
}
}

void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array<double>& Ag, const Array<double>& Yg, const Array<double>& Dg, const bool assembly)
{
using namespace consts;

// Get dimensions
const int eNoN = lM.eNoN;
const int dof = com_mod.dof;

// Get properties
auto cDmn = com_mod.cDmn;
const int cEq = com_mod.cEq;
const auto& eq = com_mod.eq[cEq];

// Initialize residual and tangent
Vector<int> ptr(eNoN);
Array<double> lR(dof, eNoN);
Array3<double> lK(dof * dof, eNoN, eNoN);

// Loop over all elements of mesh
for (int e = 0; e < lM.nEl; e++) {
// Reset
ptr = 0;
lR = 0.0;
lK = 0.0;

// Update domain and proceed if domain phys and eqn phys match
cDmn = all_fun::domain(com_mod, lM, cEq, e);
auto cPhys = eq.dmn[cDmn].phys;
if (cPhys != EquationType::phys_struct) {
continue;
}

// Update shape functions for NURBS
if (lM.eType == ElementType::NRB) {
//CALL NRBNNX(lM, e)
}

// Evaluate solid equations
eval_dsolid(e, com_mod, cep_mod, lM, Ag, Yg, Dg, ptr, lR, lK);

// Assemble into global residual and tangent
if (assembly) {
#ifdef WITH_TRILINOS
if (eq.assmTLS) {
trilinos_doassem_(const_cast<int&>(eNoN), ptr.data(), lK.data(), lR.data());
} else {
#endif
lhsa_ns::do_assem(com_mod, eNoN, ptr, lK, lR);
#ifdef WITH_TRILINOS
}
#endif
}
}
}

/// @brief Replicates the Fortan 'CONSTRUCT_dSOLID' subroutine.
//
void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array<double>& Ag, const Array<double>& Yg, const Array<double>& Dg)
void eval_dsolid(const int &e, ComMod &com_mod, CepMod &cep_mod,
const mshType &lM, const Array<double> &Ag,
const Array<double> &Yg, const Array<double> &Dg,
Vector<int> &ptr, Array<double> &lR, Array3<double> &lK)
{
using namespace consts;

Expand All @@ -216,9 +274,6 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const
const int nsd = com_mod.nsd;
const int tDof = com_mod.tDof;
const int dof = com_mod.dof;
const int cEq = com_mod.cEq;
const auto& eq = com_mod.eq[cEq];
auto cDmn = com_mod.cDmn;
const int nsymd = com_mod.nsymd;
auto& pS0 = com_mod.pS0;
auto& pSn = com_mod.pSn;
Expand All @@ -240,142 +295,106 @@ void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const
#endif

// STRUCT: dof = nsd

Vector<int> ptr(eNoN);
Vector<double> pSl(nsymd), ya_l(eNoN), N(eNoN), gr_int_g(com_mod.nGrInt), gr_props_g(lM.n_gr_props);
Array<double> xl(nsd,eNoN), al(tDof,eNoN), yl(tDof,eNoN), dl(tDof,eNoN),
bfl(nsd,eNoN), fN(nsd,nFn), pS0l(nsymd,eNoN), Nx(nsd,eNoN), lR(dof,eNoN),
bfl(nsd,eNoN), fN(nsd,nFn), pS0l(nsymd,eNoN), Nx(nsd,eNoN),
gr_props_l(lM.n_gr_props,eNoN);
Array3<double> lK(dof*dof,eNoN,eNoN);

// Create local copies
fN = 0.0;
pS0l = 0.0;
ya_l = 0.0;
gr_int_g = 0.0;
gr_props_l = 0.0;
gr_props_g = 0.0;

// Loop over all elements of mesh
for (int a = 0; a < eNoN; a++) {
int Ac = lM.IEN(a,e);
ptr(a) = Ac;

for (int e = 0; e < lM.nEl; e++) {
// Update domain and proceed if domain phys and eqn phys match
cDmn = all_fun::domain(com_mod, lM, cEq, e);
auto cPhys = eq.dmn[cDmn].phys;
if (cPhys != EquationType::phys_struct) {
continue;
for (int i = 0; i < nsd; i++) {
xl(i,a) = com_mod.x(i,Ac);
bfl(i,a) = com_mod.Bf(i,Ac);
}

// Update shape functions for NURBS
if (lM.eType == ElementType::NRB) {
//CALL NRBNNX(lM, e)
for (int i = 0; i < tDof; i++) {
al(i,a) = Ag(i,Ac);
dl(i,a) = Dg(i,Ac);
yl(i,a) = Yg(i,Ac);
}

// Create local copies
fN = 0.0;
pS0l = 0.0;
ya_l = 0.0;
gr_int_g = 0.0;
gr_props_l = 0.0;
gr_props_g = 0.0;

for (int a = 0; a < eNoN; a++) {
int Ac = lM.IEN(a,e);
ptr(a) = Ac;

for (int i = 0; i < nsd; i++) {
xl(i,a) = com_mod.x(i,Ac);
bfl(i,a) = com_mod.Bf(i,Ac);
}

for (int i = 0; i < tDof; i++) {
al(i,a) = Ag(i,Ac);
dl(i,a) = Dg(i,Ac);
yl(i,a) = Yg(i,Ac);
}

if (lM.fN.size() != 0) {
for (int iFn = 0; iFn < nFn; iFn++) {
for (int i = 0; i < nsd; i++) {
fN(i,iFn) = lM.fN(i+nsd*iFn,e);
}
if (lM.fN.size() != 0) {
for (int iFn = 0; iFn < nFn; iFn++) {
for (int i = 0; i < nsd; i++) {
fN(i,iFn) = lM.fN(i+nsd*iFn,e);
}
}
}

if (pS0.size() != 0) {
pS0l.set_col(a, pS0.col(Ac));
}
if (pS0.size() != 0) {
pS0l.set_col(a, pS0.col(Ac));
}

if (cem.cpld) {
ya_l(a) = cem.Ya(Ac);
}
if (cem.cpld) {
ya_l(a) = cem.Ya(Ac);
}

if (lM.gr_props.size() != 0) {
for (int igr = 0; igr < lM.n_gr_props; igr++) {
gr_props_l(igr,a) = lM.gr_props(igr,Ac);
}
if (lM.gr_props.size() != 0) {
for (int igr = 0; igr < lM.n_gr_props; igr++) {
gr_props_l(igr,a) = lM.gr_props(igr,Ac);
}
}
}

// Gauss integration
//
lR = 0.0;
lK = 0.0;

double Jac{0.0};
Array<double> ksix(nsd,nsd);
// Gauss integration
double Jac{0.0};
Array<double> ksix(nsd,nsd);

for (int g = 0; g < lM.nG; g++) {
if (g == 0 || !lM.lShpF) {
auto Nx_g = lM.Nx.slice(g);
nn::gnn(eNoN, nsd, nsd, Nx_g, xl, Nx, Jac, ksix);
if (utils::is_zero(Jac)) {
throw std::runtime_error("[construct_dsolid] Jacobian for element " + std::to_string(e) + " is < 0.");
}
for (int g = 0; g < lM.nG; g++) {
if (g == 0 || !lM.lShpF) {
auto Nx_g = lM.Nx.slice(g);
nn::gnn(eNoN, nsd, nsd, Nx_g, xl, Nx, Jac, ksix);
if (utils::is_zero(Jac)) {
throw std::runtime_error("[construct_dsolid] Jacobian for element " + std::to_string(e) + " is < 0.");
}
double w = lM.w(g) * Jac;
N = lM.N.col(g);
pSl = 0.0;

// Get internal growth and remodeling variables
if (com_mod.grEq) {
// todo mrp089: add a function like rslice for vectors to Array3
for (int i = 0; i < com_mod.nGrInt; i++) {
gr_int_g(i) = com_mod.grInt(e,g,i);
}
}
double w = lM.w(g) * Jac;
N = lM.N.col(g);
pSl = 0.0;

// Get internal growth and remodeling variables
if (com_mod.grEq) {
// todo mrp089: add a function like rslice for vectors to Array3
for (int i = 0; i < com_mod.nGrInt; i++) {
gr_int_g(i) = com_mod.grInt(e,g,i);
}
}

if (nsd == 3) {
struct_3d(com_mod, cep_mod, eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, pS0l, pSl, ya_l, gr_int_g, gr_props_l, lR, lK);
if (nsd == 3) {
struct_3d(com_mod, cep_mod, eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, pS0l, pSl, ya_l, gr_int_g, gr_props_l, lR, lK);

} else if (nsd == 2) {
struct_2d(com_mod, cep_mod, eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, pS0l, pSl, ya_l, gr_int_g, gr_props_l, lR, lK);
}
} else if (nsd == 2) {
struct_2d(com_mod, cep_mod, eNoN, nFn, w, N, Nx, al, yl, dl, bfl, fN, pS0l, pSl, ya_l, gr_int_g, gr_props_l, lR, lK);
}

// Set internal growth and remodeling variables
if (com_mod.grEq) {
// todo mrp089: add a function like rslice for vectors to Array3
for (int i = 0; i < com_mod.nGrInt; i++) {
com_mod.grInt(e,g,i) = gr_int_g(i);
}
// Set internal growth and remodeling variables
if (com_mod.grEq) {
// todo mrp089: add a function like rslice for vectors to Array3
for (int i = 0; i < com_mod.nGrInt; i++) {
com_mod.grInt(e,g,i) = gr_int_g(i);
}
}

// Prestress
if (pstEq) {
for (int a = 0; a < eNoN; a++) {
int Ac = ptr(a);
pSa(Ac) = pSa(Ac) + w*N(a);
for (int i = 0; i < pSn.nrows(); i++) {
pSn(i,Ac) = pSn(i,Ac) + w*N(a)*pSl(i);
}
// Prestress
if (pstEq) {
for (int a = 0; a < eNoN; a++) {
int Ac = ptr(a);
pSa(Ac) = pSa(Ac) + w*N(a);
for (int i = 0; i < pSn.nrows(); i++) {
pSn(i,Ac) = pSn(i,Ac) + w*N(a)*pSl(i);
}
}
}

// Assembly
//
#ifdef WITH_TRILINOS
if (eq.assmTLS) {
trilinos_doassem_(const_cast<int&>(eNoN), ptr.data(), lK.data(), lR.data());
} else {
#endif
lhsa_ns::do_assem(com_mod, eNoN, ptr, lK, lR);
#ifdef WITH_TRILINOS
}
#endif

}
}

Expand Down
5 changes: 4 additions & 1 deletion Code/Source/svFSI/sv_struct.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,10 @@ void b_struct_3d(const ComMod& com_mod, const int eNoN, const double w, const Ve
Array<double>& lR, Array3<double>& lK);

void construct_dsolid(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array<double>& Ag,
const Array<double>& Yg, const Array<double>& Dg);
const Array<double>& Yg, const Array<double>& Dg, const bool assembly=true);

void eval_dsolid(const int& e, ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Array<double>& Ag,
const Array<double>& Yg, const Array<double>& Dg, Vector<int>& ptr, Array<double>& lR, Array3<double>& lK);

void struct_2d(ComMod& com_mod, CepMod& cep_mod, const int eNoN, const int nFn, const double w,
const Vector<double>& N, const Array<double>& Nx, const Array<double>& al, const Array<double>& yl,
Expand Down

0 comments on commit a445e1f

Please sign in to comment.