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/liner_solver/add_bc_mul.cpp

@@ -55,6 +55,9 @@ void add_bc_mul(FSILS_lhsType& lhs, const BcopType op_Type, const int dof, const
   Vector<double> coef(lhs.nFaces); 
   Array<double> v(dof,lhs.nNo);
 
+  //Setting coef depending on adding resistance to stiffness or
+  //computing preconditioner
+
   if (op_Type == BcopType::BCOP_TYPE_ADD) {
     for (int i = 0; i < lhs.nFaces; i++) {
       coef(i) = lhs.face[i].res;
@@ -70,6 +73,25 @@ void add_bc_mul(FSILS_lhsType& lhs, const BcopType op_Type, const int dof, const
 
   for (int faIn = 0; faIn < lhs.nFaces; faIn++) {
     auto& face = lhs.face[faIn];
+    // Cap faces do not contribute to the tangent matrix
+    if (face.isCap) {
+      continue;
+    }
+
+    // In the following calculations, we are computing the product of the
+    // coupled BC tangent contribution with the vector X (refer to Moghadam
+    // et al. 2013 eq. 27). This is computed by first constructing the vector
+    // v, which one of the integrals found in the expression, int{N_A * n_i} dGamma.
+    // Then, v is dotted with X to yield a quantity S. Then S is multiplied by
+    // by v again, and also multiplied by the appropriate coefficients in
+    // the expression.
+    // The calculations are complicated somewhat if there is a capping surface,
+    // but these complications are explained below.
+
+
+    // Calculating S, which is the inner product of the right integral (v) and
+    // the vector to be multiplied (X).
+
     int nsd = std::min(face.dof, dof);
 
     if (face.coupledFlag) {
@@ -85,15 +107,32 @@ void add_bc_mul(FSILS_lhsType& lhs, const BcopType op_Type, const int dof, const
         }
         // Computing S = coef * v^T * X
         double S = coef(faIn) * dot::fsils_dot_v(dof, lhs.mynNo, lhs.commu, v, X);
+
+        // Add cap surface contribution to S
+        // Cap surfaces contribute to flow rate but not pressure
+        if (face.isCapped) {
+          int faInCap = face.faInCap;
+          auto& faceCap = lhs.face[faInCap];
 
+          Array<double> vcap(dof,lhs.nNo);
+          vcap = 0.0;
+          for (int a = 0; a < faceCap.nNo; a++) {
+            int Ac = faceCap.glob(a);
+            for (int i = 0; i < nsd; i++) {
+              vcap(i,Ac) = faceCap.valM(i,a);  
+            }
+          }
+          // Add capping surface contribution to S
+          S = S + coef(faIn) * dot::fsils_dot_v(dof, lhs.mynNo, lhs.commu, vcap, X);
+        }
+
         // Computing Y = Y + v * S
         for (int a = 0; a < face.nNo; a++) {
           int Ac = face.glob(a);
           for (int i = 0; i < nsd; i++) {
             Y(i,Ac) = Y(i,Ac) + v(i,Ac)*S;
           }
         }
-
       } 
       // If face is not shared across procs
       else  {
@@ -105,8 +144,29 @@ void add_bc_mul(FSILS_lhsType& lhs, const BcopType op_Type, const int dof, const
             S = S + face.valM(i,a)*X(i,Ac);
           }
         }
-        S = coef(faIn) * S;
 
+        // If capping surface is present add its contribution to S
+        if (face.isCapped) {
+          int faInCap = face.faInCap;
+          auto& faceCap = lhs.face[faInCap];
+
+          if (!faceCap.coupledFlag) {
+            std::cerr << "ADDBCMUL(): Cap face is not coupled. Probably cap face has zero resistance." << std::endl;
+            throw std::runtime_error("FSILS: FATAL ERROR");
+          }
+
+          for (int a = 0; a < faceCap.nNo; a++) {
+            int Ac = faceCap.glob(a);
+            for (int i = 0; i < nsd; i++) {
+              S = S + faceCap.valM(i,a)*X(i,Ac);
+            }
+          }
+        }
+
+        // Multiply S by the resistance or related quantity if
+        // preconditioning
+        S = coef(faIn) * S;
+
         // Computing Y = Y + v * S
         for (int a = 0; a < face.nNo; a++) {
           int Ac = face.glob(a);

Code/Source/liner_solver/bc.cpp

@@ -43,7 +43,7 @@ namespace fsi_linear_solver {
 ///  lhs.face[faIn].valM
 //
 void fsils_bc_create(FSILS_lhsType& lhs, int faIn, int nNo, int dof, BcType BC_type, const Vector<int>& gNodes, 
-    const Array<double>& Val)
+    const Array<double>& Val, bool isCap)
 {
   using namespace consts;
 
@@ -60,6 +60,7 @@ void fsils_bc_create(FSILS_lhsType& lhs, int faIn, int nNo, int dof, BcType BC_t
   dmsg << "Val.size(): " << Val.size();
   dmsg << "Val.nrows: " << Val.nrows_;
   dmsg << "Val.ncols: " << Val.ncols_;
+  dmsg << "isCap: " << isCap;
   #endif
 
   if (faIn >= lhs.nFaces) {
@@ -71,9 +72,12 @@ void fsils_bc_create(FSILS_lhsType& lhs, int faIn, int nNo, int dof, BcType BC_t
     throw std::runtime_error("FSILS: faIn is smaller than zero");
   }
 
+  lhs.face[faIn].foC = true;
   lhs.face[faIn].nNo = nNo;
   lhs.face[faIn].dof = dof;
   lhs.face[faIn].bGrp = BC_type;
+  // Set cap flag for face
+  lhs.face[faIn].isCap = isCap;
 
   lhs.face[faIn].glob.resize(nNo); 
   lhs.face[faIn].val.resize(dof,nNo);   
@@ -116,7 +120,6 @@ void fsils_bc_create(FSILS_lhsType& lhs, int faIn, int nNo, int dof, BcType BC_t
           v(i,Ac) = lhs.face[faIn].val(i,a);
         }
       }
-
       fsils_commuv(lhs, dof, v); 
 
       for (int a = 0; a < nNo; a++) {
@@ -184,9 +187,8 @@ void fsils_bc_update(FSILS_lhsType& lhs, int faIn, int nNo, int dof, const Array
   }
 
   // Communicate update among procs
-  if (lhs.face[faIn].sharedFlag){
+  if (lhs.face[faIn].sharedFlag && !lhs.face[faIn].isCap) {
     Array<double> v(dof,lhs.nNo);
-    v = 0.0;
 
     for (int a = 0; a < nNo; a++) {
       int Ac = lhs.face[faIn].glob(a);

Code/Source/liner_solver/fils_struct.hpp

@@ -154,6 +154,21 @@ class FSILS_faceType
 
     /// Neu W*Sai                      (TMP)
     Array<double> valM;
+
+    /// Flag for cap face          (USE)
+    bool isCap = false;
+
+    /// Index of mesh in msh associated with this face
+    int iM = -1;
+
+    /// Index of face in msh[iM].fa that caps this face
+    int iFa = -1;
+
+    /// Index of face in lhs.face that caps this face
+    int faInCap = -1;
+
+    /// Flag indicating if this face is capped by another face
+    bool isCapped = false;
 };
 
 /// @brief Modified in:

Code/Source/liner_solver/fsils_api.hpp

@@ -39,7 +39,7 @@
 namespace fsi_linear_solver {
 
 void fsils_bc_create(FSILS_lhsType& lhs, int faIn, int nNo, int dof, BcType BC_type, const Vector<int>& gNodes,
-    const Array<double>& Val);
+    const Array<double>& Val, bool isCap = false);
 
 void fsils_bc_create(FSILS_lhsType& lhs, int faIn, int nNo, int dof, BcType BC_type, const Vector<int>& gNodes);
 

Code/Source/solver/ComMod.h

@@ -73,7 +73,7 @@ class fcType
 {
   public:
 
-    bool defined() { return n != 0; };
+    bool defined() const { return n != 0; };
 
     // If this is a ramp function
     bool lrmp = false;
@@ -109,7 +109,7 @@ class MBType
 {
   public:
 
-    bool defined() { return dof != 0; };
+    bool defined() const { return dof != 0; };
 
     // Degrees of freedom of d(:,.,.)
     int dof = 0;
@@ -180,6 +180,15 @@ class bcType
     // Pointer to FSILS%bc
     int lsPtr = -1;
 
+    // Index of cap BC associated wtih this BC
+    int iCapBC = -1;
+
+    // Flag indicating if this BC has a capping BC
+    bool hasCapBC = false;
+
+    // Name of face that caps this surface
+    std::string capName;
+
     // Undeforming Neu BC master-slave node parameters.
     int masN = 0;
 
@@ -544,6 +553,9 @@ class faceType
 
     //faceType& operator=(const faceType& rhs);
 
+    // Flag for cap face (i.e. face does not lie on volume mesh)
+    bool isCap = false;
+
     // Parametric direction normal to this face (NURBS)
     int d = 0;
 
@@ -571,6 +583,9 @@ class faceType
     // Number of nodes
     int nNo = 0;
 
+    //ID number of cap face that caps this face
+    int capID = -1;
+
     // Global element Ids
     Vector<int> gE;
 
@@ -919,6 +934,9 @@ class mshType
     /// @brief Whether the mesh is fibers (Purkinje)
     bool lFib = false;
 
+    /// @brief Whether the mesh is a cap
+    bool isCap = false;
+
     /// @brief Element type
     consts::ElementType eType = consts::ElementType::NA;
     //int eType = eType_NA

Code/Source/solver/Parameters.cpp

@@ -459,6 +459,9 @@ BoundaryConditionParameters::BoundaryConditionParameters()
   set_parameter("Zero_out_perimeter", false, !required, zero_out_perimeter);
 
   set_parameter("Resistance", 1.e5, !required, resistance);
+
+  set_parameter("Capping_face", "", !required, capping_face);
+  set_parameter("cap_svZeroDSolver_block", "", !required, cap_svzerod_solver_block);
 }
 
 void BoundaryConditionParameters::print_parameters()
@@ -2275,6 +2278,7 @@ FaceParameters::FaceParameters()
 
   set_parameter("End_nodes_face_file_path", "", !required, end_nodes_face_file_path);
   set_parameter("Face_file_path", "", !required, face_file_path);
+  set_parameter("Cap", false, !required, cap_face);
 
   set_parameter("Quadrature_modifier_TRI3", (2.0/3.0), !required, quadrature_modifier_TRI3);
 }

Code/Source/solver/Parameters.h

@@ -821,6 +821,9 @@ class BoundaryConditionParameters : public ParameterLists
     Parameter<bool> zero_out_perimeter;
 
     Parameter<double> resistance;
+
+    Parameter<std::string> capping_face;
+    Parameter<std::string> cap_svzerod_solver_block;
 };
 
 /// @brief The OutputParameters class stores parameters for the
@@ -1462,6 +1465,7 @@ class FaceParameters : public ParameterLists
     Parameter<std::string> end_nodes_face_file_path;
     Parameter<std::string> face_file_path;
     Parameter<std::string> name;
+    Parameter<bool> cap_face;
 
     Parameter<double> quadrature_modifier_TRI3;
 };