3 d0 d coupling cap

.gitignore

@@ -35,3 +35,6 @@ Q_svZeroD
 P_svZeroD
 # local svZeroD build
 svZeroDSolver
+
+# test artifacts (VTK, processed results)
+tests/**/calc_volume_struct_output/

Code/Source/solver/Parameters.cpp

@@ -459,6 +459,8 @@ BoundaryConditionParameters::BoundaryConditionParameters()
   set_parameter("Weakly_applied", false, !required, weakly_applied);
   set_parameter("Zero_out_perimeter", false, !required, zero_out_perimeter);
 
+  set_parameter("Resistance", 1.e5, !required, resistance);
+
   set_parameter("Capping_face", "", !required, capping_face);
 }
 

Code/Source/solver/Parameters.h

@@ -821,6 +821,8 @@ class BoundaryConditionParameters : public ParameterLists
     Parameter<bool> weakly_applied;
     Parameter<bool> zero_out_perimeter;
 
+    Parameter<double> resistance;
+
     Parameter<std::string> capping_face;
 };
 

Code/Source/solver/all_fun.cpp

@@ -437,6 +437,153 @@ global(const ComMod& com_mod, const CmMod& cm_mod, const mshType& lM, const Arra
   return result;
 }
 
+///////////////////////////////////////////////////////////////////////////
+/// Convert the vIntegM function from svFSI to c++.
+/// @brief This routine integrated a scalar field over a particular domain.
+///
+/// @param iM mesh index
+/// @param s an array containing a scalar value for each node in the mesh
+/// Replicates 'FUNCTION vIntegM(dId, s, l, u, pFlag)' defined in ALLFUN.f.
+//
+double integ(const ComMod& com_mod, const CmMod& cm_mod, int iM, const Array<double>& s)
+{
+  using namespace consts;
+
+  #define n_debug_integ_v
+  #ifdef debug_integ_v
+  DebugMsg dmsg(__func__, com_mod.cm.idcm());
+  dmsg.banner();
+  dmsg << "vIntegM " << " ";
+  dmsg << "iM: " << iM;
+  #endif
+
+  int nNo = s.ncols();
+  int tnNo = com_mod.tnNo;
+  bool ibFlag = com_mod.ibFlag;
+
+  if (nNo != tnNo) {
+    std::string msg = "Incompatible vector size in vIntegM in domain: ";
+    msg += std::to_string(iM);
+    throw std::runtime_error(msg);
+  }
+
+  double result = 0.0; 
+  int nsd = com_mod.nsd;
+  fsType fs;
+
+  auto& msh = com_mod.msh[iM];
+  int insd = nsd;
+  if (msh.lShl) insd = nsd-1;
+  if (msh.lFib) insd = 1;
+
+  // Update pressure function space for Taylor-Hood type element
+  if (com_mod.msh[iM].nFs == 2) {
+    fs.nG    = msh.fs[1].nG;
+    fs.eType = msh.fs[1].eType;
+    fs.lShpF = msh.fs[1].lShpF;
+    fs.eNoN  = msh.fs[1].eNoN;
+
+    fs.w.resize(fs.nG); 
+    fs.N.resize(fs.eNoN,fs.nG); 
+    fs.Nx.resize(nsd,fs.eNoN,fs.nG);
+
+    if (fs.eType != ElementType::NRB) {
+      fs.w  = msh.fs[1].w;
+      fs.N  = msh.fs[1].N;
+      fs.Nx = msh.fs[1].Nx;
+    }
+
+  } else { 
+    fs.nG    = msh.fs[0].nG;
+    fs.eType = msh.fs[0].eType;
+    fs.lShpF = msh.fs[0].lShpF;
+    fs.eNoN  = msh.fs[0].eNoN;
+
+    fs.w.resize(fs.nG); 
+    fs.N.resize(fs.eNoN,fs.nG); 
+    fs.Nx.resize(nsd,fs.eNoN,fs.nG);
+
+    if (fs.eType != ElementType::NRB) {
+      fs.w  = msh.fs[0].w;
+      fs.N  = msh.fs[0].N;
+      fs.Nx = msh.fs[0].Nx;
+    }
+  }
+  int eNoN = fs.eNoN;
+
+  Array<double> xl(nsd,eNoN); 
+  Array<double> Nxi(insd,eNoN); 
+  Array<double> Nx(insd,eNoN); 
+  Vector<double> sl(eNoN); 
+  Array<double> tmps(nsd,insd);
+  Array<double> tmp(nsd,nsd);
+
+  for (int e = 0; e < msh.nEl; e++) {
+    // Updating the shape functions, if this is a NURB
+    //
+    // [TODO:DaveP] not implemented. 
+    //
+    if (msh.eType == ElementType::NRB) {
+      //CALL NRBNNX(msh(iM), e)
+      fs.w  = msh.w;
+      fs.N  = msh.N;
+      fs.Nx = msh.Nx;
+    }
+
+    int ibl = 0;
+    for (int a = 0; a < eNoN; a++) { 
+      int Ac = msh.IEN(a,e);
+      xl.set_col(a, com_mod.x.col(Ac));
+
+      if (com_mod.mvMsh) {
+        for (int i = 0; i < nsd; i++) { 
+          xl(i,a) += com_mod.Do(i+nsd+1,Ac);
+        }
+      }
+      sl(a) = s(0,Ac);
+      ibl = ibl + com_mod.iblank(Ac);
+    }
+
+    if (ibl == eNoN) {
+      continue;
+    }
+
+    double Jac = 0.0;
+
+    for (int g = 0; g < fs.nG; g++) {
+      Nxi = fs.Nx.slice(g);
+
+      if (g == 0 || !fs.lShpF) {
+        if (msh.lShl) {
+          Vector<double> nV(nsd);
+          nn::gnns(nsd, eNoN, Nxi, xl, nV, tmps, tmps);
+          Jac  = sqrt(utils::norm(nV));
+        } else { 
+          nn::gnn(eNoN, nsd, insd, Nxi, xl, Nx, Jac, tmp);
+        }
+      }
+
+      if (utils::is_zero(Jac)) {
+        throw std::runtime_error("Jac < 0 for element: " + std::to_string(e) + ")");
+      }
+      double sHat = 0.0;
+
+      for (int a = 0; a < eNoN; a++) {
+        int Ac = msh.IEN(a,e);
+        sHat = sHat + sl(a)*fs.N(a,g);
+      }
+      result += fs.w(g) * Jac * sHat;
+    }
+  }
+
+  if (com_mod.cm.seq()) {
+    return result;
+  }
+  result = com_mod.cm.reduce(cm_mod, result);
+  return result;
+}
+///////////////////////////////////////////////////////////////////////////
+
 /// @brief This routine integrates a scalar field over a particular domain.
 ///
 /// Note that 'l' and 'u' should be 0-based and are used to index into 's'.

Code/Source/solver/set_bc.cpp

@@ -1370,11 +1370,8 @@ void set_bc_neu(ComMod& com_mod, const CmMod& cm_mod, const Array<double>& Yg, c
     dmsg << "----- iBc " << iBc+1;
     #endif
 
-    // Skip if the face is virtual
-    if (com_mod.msh[iM].fa[iFa].vrtual && utils::btest(bc.bType, iBC_res) && bc.flwP) {
-      eq_assem::fsi_ls_upd(com_mod, cm_mod, bc, com_mod.msh[iM].fa[iFa]);
-      continue;
-    }
+    if (utils::btest(bc.bType, iBC_Ris0D))  {continue;}
+
     if (utils::btest(bc.bType, iBC_Neu)) {
       #ifdef debug_set_bc_neu
       dmsg << "iM: " << iM+1;
@@ -1499,15 +1496,14 @@ void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const
     }
   }
   // Now treat Robin BC (stiffness and damping) here
-  //
-  if (utils::btest(lBc.bType,iBC_Robin)) {
-    set_bc_rbnl(com_mod, cm_mod, lFa, lBc.k, lBc.c, lBc.rbnN, Yg, Dg);
+  if (lBc.robin_bc.is_initialized()) {
+    set_bc_rbnl(com_mod, cm_mod, lFa, lBc.robin_bc, Yg, Dg);
   }
 }
 
 /// @brief Set Robin BC contribution to residual and tangent
 //
-void set_bc_rbnl(ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const double ks, const double cs, const bool isN, 
+void set_bc_rbnl(ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const RobinBoundaryCondition& robin_bc,
   const Array<double>& Yg, const Array<double>& Dg)
 {
   using namespace consts;

Code/Source/solver/set_bc.h

@@ -62,7 +62,7 @@ void set_bc_dir_wl(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, cons
 void set_bc_neu(ComMod& com_mod, const CmMod& cm_mod, const Array<double>& Yg, const Array<double>& Dg);
 void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const faceType& lFa, const Array<double>& Yg, const Array<double>& Dg);
 
-void set_bc_rbnl(ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const double ks, const double cs, const bool isN, 
+void set_bc_rbnl(ComMod& com_mod, const CmMod& cm_mod, const faceType& lFa, const RobinBoundaryCondition& robin_bc,
   const Array<double>& Yg, const Array<double>& Dg);
 
 void set_bc_trac_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const faceType& lFa);

tests/cases/struct/LV_NeoHookean_passive_genBC/process_results.py

@@ -43,7 +43,7 @@
 
 # svFSI results folder, containing results.vtu
 #results_folder = os.path.join(sim_folder, 'results_svfsi/')
-results_folder = os.path.join(sim_folder, '4-procs/')
+results_folder = os.path.join(sim_folder, '6-procs/')
 
 # File containing genBC output
 #alldata_file = 'AllData_svfsi_vvedula22'