Initial pass at using linearly interpolated areas within the solver using the spatial characteristics. This is failing system tests so clearly there are issues. I have a few notes on things to check but mostly I just need to go back and verify all the locations where we're using the new data. If that isn't sufficient I might need to do some spot checks by outputing what the solver is seeing for those values.

Author: arsLibera

Code/Source/cvOneDBFSolver.cxx

@@ -95,6 +95,15 @@ double                        cvOneDBFSolver::convCriteria = 0;
 BoundCondType                 cvOneDBFSolver::inletBCtype;
 int                           cvOneDBFSolver::ASCII = 1;
 
+
+namespace{
+
+  double linearInterpolate(double x, double x1, double y1, double x2, double y2) {
+      return y1 + (x - x1) * (y2 - y1) / (x2 - x1);
+  }
+      
+} // namespace
+
 // SET MODE PTR
 void cvOneDBFSolver::SetModelPtr(cvOneDModel *mdl){
   model = mdl;
@@ -424,7 +433,6 @@ void cvOneDBFSolver::postprocess_VTK_XML3D_ONEFILE(){
   cvDoubleVec tmp;
   cvDoubleMat segNodeList;
   double lengthByNodes = 0.0;
-  double lengthBySegment = 0.0;
   int startOut = 0;
   int finishOut = 0;
   double segLength = 0.0;
@@ -470,7 +478,6 @@ void cvOneDBFSolver::postprocess_VTK_XML3D_ONEFILE(){
     segVers[2][0] /= mod;
 
     lengthByNodes = mod;
-    lengthBySegment = currSeg->getSegmentLength();
 
     // Compute Segment Local axis system
     evalSegmentLocalAxis(segVers);
@@ -484,8 +491,12 @@ void cvOneDBFSolver::postprocess_VTK_XML3D_ONEFILE(){
       currCentre[2] = nodeList[inletSegJoint][2] + loopEl*lengthByNodes/double(currSeg->getNumElements())*segVers[2][0];
 
       // Get initial radius at current location
-      currIniArea = currSeg->getInitInletS() + (loopEl/double(currSeg->getNumElements()))*(currSeg->getInitOutletS() - currSeg->getInitInletS());
-      currIniRad = sqrt(currIniArea/M_PI);
+
+      // TODO: verify this is the correct "z" location we want
+      // to interpolate at. These loops and myriad variables make
+      // something super simple super confusing.
+      double const zAxial = linearInterpolate(loopEl, 0, currSeg->getInletZ(), currSeg->getNumElements(), currSeg->getOutletZ());
+      currIniRad = currSeg->getInitialRadius(zAxial);
 
       // Loop on the subdivisions
       for(int loopSubdiv=0;loopSubdiv<circSubdiv;loopSubdiv++){
@@ -559,7 +570,14 @@ void cvOneDBFSolver::postprocess_VTK_XML3D_ONEFILE(){
       for(int j=startOut;j<finishOut;j+=2){
 
         // Evaluate Initial Area at current location
-        iniArea = currSeg->getInitInletS() + (((j-startOut)/2)/double(currSeg->getNumElements()))*(currSeg->getInitOutletS() - currSeg->getInitInletS());
+        
+        // TODO: verify this is the correct "z" location we want
+        // to interpolate at. These loops and myriad variables make
+        // something super simple super confusing.
+        double const zAxial = linearInterpolate((j-startOut)/2.0, 
+          startOut, currSeg->getInletZ(), finishOut, currSeg->getOutletZ());
+        iniArea = currSeg->getInitialArea(zAxial);
+
         // Eval Current Area at current location
         newArea = TotalSolution[loopTime][j];
         // Evaluate Radial displacement
@@ -738,7 +756,6 @@ void cvOneDBFSolver::postprocess_VTK_XML3D_MULTIPLEFILES(){
   cvDoubleVec tmp;
   cvDoubleMat segNodeList;
   double lengthByNodes = 0.0;
-  double lengthBySegment = 0.0;
   int startOut = 0;
   int finishOut = 0;
   double segLength = 0.0;
@@ -812,7 +829,6 @@ void cvOneDBFSolver::postprocess_VTK_XML3D_MULTIPLEFILES(){
       segVers[2][0] /= mod;
 
       lengthByNodes = mod;
-      lengthBySegment = currSeg->getSegmentLength();
 
       // Compute Segment Local axis system
       evalSegmentLocalAxis(segVers);
@@ -826,9 +842,13 @@ void cvOneDBFSolver::postprocess_VTK_XML3D_MULTIPLEFILES(){
         currCentre[2] = nodeList[inletSegJoint][2] + loopEl*lengthByNodes/double(currSeg->getNumElements())*segVers[2][0];
 
         // Get initial radius at current location
-        currIniArea = currSeg->getInitInletS() + (loopEl/double(currSeg->getNumElements()))*(currSeg->getInitOutletS() - currSeg->getInitInletS());
-        currIniRad = sqrt(currIniArea/M_PI);
-
+             
+        // TODO: verify this is the correct "z" location we want
+        // to interpolate at. These loops and myriad variables make
+        // something super simple super confusing.
+        double const zAxial = linearInterpolate(loopEl, 0, currSeg->getInletZ(), currSeg->getNumElements(), currSeg->getOutletZ());
+        currIniRad = currSeg->getInitialRadius(zAxial);
+        
         // Loop on the subdivisions
         for(int loopSubdiv=0;loopSubdiv<circSubdiv;loopSubdiv++){
           currTheta = loopSubdiv*2*M_PI/double(circSubdiv);
@@ -899,7 +919,14 @@ void cvOneDBFSolver::postprocess_VTK_XML3D_MULTIPLEFILES(){
       for(int j=startOut;j<finishOut;j+=2){
 
         // Evaluate Initial Area at current location
-        iniArea = currSeg->getInitInletS() + (((j-startOut)/2)/double(currSeg->getNumElements()))*(currSeg->getInitOutletS() - currSeg->getInitInletS());
+             
+        // TODO: verify this is the correct "z" location we want
+        // to interpolate at. These loops and myriad variables make
+        // something super simple super confusing.
+        double const zAxial = linearInterpolate((j-startOut)/2.0, 
+          startOut, currSeg->getInletZ(), finishOut, currSeg->getOutletZ());
+        iniArea = currSeg->getInitialArea(zAxial);
+        
         // Eval Current Area at current location
         newArea = TotalSolution[loopTime][j];
         // Evaluate Radial displacement
@@ -1124,19 +1151,15 @@ void cvOneDBFSolver::QuerryModelInformation(void)
     for (i=0; i<is; i++){
       cvOneDSegment* seg = model->getSegment(i);
       long nels = seg->getNumElements();
-      double segLen = seg->getSegmentLength();
       int matID = seg->getMaterialID();
       MeshType mType = seg->getMeshType();
-      double zin = seg->getInletZ();
-      double zout = seg->getOutletZ();
 
       cvOneDSubdomain* subdomain = new cvOneDSubdomain;
       assert(subdomain != 0);
       subdomain -> SetNumberOfNodes(nels+1);
       subdomain -> SetNumberOfElements(nels);
       subdomain -> SetMeshType(mType);
-      subdomain -> Init(zin, zout);
-
+      subdomain -> Init(seg->getSpatialCharacteristics());
 
       // Get the Initial Properties of the subdomain...
       double Qo = 0.0;
@@ -1147,17 +1170,13 @@ void cvOneDBFSolver::QuerryModelInformation(void)
         P0 =  seg->getInitialPressure();
         dQ0_dT = 0.0;
       }
-      double So = seg->getInitInletS();
-      double Sn = seg->getInitOutletS();
       BoundCondType boundT = seg -> getBoundCondition();
       double  boundV= seg -> getBoundValue();
 
       // Set these in the subdomain.
       subdomain->SetInitialFlow(Qo);
       subdomain->SetInitialdFlowdT(dQ0_dT);
-      subdomain->SetInitInletS(So);
       subdomain->SetInitialPressure(P0);
-      subdomain->SetInitOutletS(Sn);
       subdomain->SetGlobal1stNodeID(temp);
       subdomain->SetBoundCondition(boundT);
       if(!seg->IsOutlet){
@@ -1324,22 +1343,21 @@ void cvOneDBFSolver::CreateGlobalArrays(void){
 
 // Initialize the solution, that is, area as area input and flow rate as 0 except the inlet
 void cvOneDBFSolver::CalcInitProps(long ID){
-  double segLen = subdomainList[ID] -> GetLength();
   double Qo, dQ0dT;
-  Qo = subdomainList[ID] -> GetInitialFlow();
+  auto const& subdomain = subdomainList[ID];
+  Qo = subdomain -> GetInitialFlow();
   dQ0dT=0;
 
-  double So = subdomainList[ID] -> GetInitInletS();
-  double Sn = subdomainList[ID] -> GetInitOutletS();
-  for( long node = 0; node < subdomainList[ID]->GetNumberOfNodes(); node++){
-  double zn = subdomainList[ID]->GetNodalCoordinate( node);
+  for( long node = 0; node < subdomain->GetNumberOfNodes(); node++){
+  double zn = subdomain->GetNodalCoordinate( node);
   long eqNumbers[2];
   mathModels[0]->GetNodalEquationNumbers(node, eqNumbers, ID);
 
   // Linear Interpolation
-  double zi = (zn - segLen)/(0.0-segLen);
-  double Si = (zi*(So - Sn)) + Sn;
-  (*previousSolution)[eqNumbers[0]] = Si;
+  // This is "Si"
+  // TODO: verify this calcluation makes sense
+  double initialAreaAtZn = getInterpolatedArea(zn, subdomain->getSpatialCharacteristics());
+  (*previousSolution)[eqNumbers[0]] = initialAreaAtZn;
 
     if(node == 0){
       (*previousSolution)[eqNumbers[1]] = Qo;

Code/Source/cvOneDMaterial.h

@@ -1,34 +1,34 @@
-/* Copyright (c) Stanford University, The Regents of the University of
- *               California, and others.
- *
- * All Rights Reserved.
- *
- * See Copyright-SimVascular.txt for additional details.
- *
- * Permission is hereby granted, free of charge, to any person obtaining
- * a copy of this software and associated documentation files (the
- * "Software"), to deal in the Software without restriction, including
- * without limitation the rights to use, copy, modify, merge, publish,
- * distribute, sublicense, and/or sell copies of the Software, and to
- * permit persons to whom the Software is furnished to do so, subject
- * to the following conditions:
- *
- * The above copyright notice and this permission notice shall be included
- * in all copies or substantial portions of the Software.
- *
- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
- * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
- * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
- * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
- * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
- * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
- * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
- * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
- * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
- * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- */
-
+/* Copyright (c) Stanford University, The Regents of the University of
+ *               California, and others.
+ *
+ * All Rights Reserved.
+ *
+ * See Copyright-SimVascular.txt for additional details.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+ * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+ * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+ * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
 #ifndef CVONEDMATERIAL_H
 #define CVONEDMATERIAL_H
 
@@ -38,8 +38,10 @@
 //  This class maintains Material Properties of the subdomain.
 //
 
-# include "cvOneDEnums.h"
-# include <math.h>
+#include "cvOneDEnums.h"
+#include "cvOneDSegmentSpatialCharacteristics.h"
+
+#include <math.h>
 
 class cvOneDMaterial{
 
@@ -82,16 +84,27 @@ class cvOneDMaterial{
     virtual double GetProperty(char* what) const = 0;
     virtual double GetIntegralpS (double area, double z) const = 0;
 
-
-    virtual double GetN(double S) const = 0;//not really dependent on S actually IV 080703
+
+    virtual double GetN(double S) const = 0;//not really dependent on S actually IV 080703
     virtual double GetEHR(double z) const = 0;
-
+
     virtual double GetOutflowFunction(double pressure, double z) const = 0;
     virtual double GetDpDz(double area, double z) const = 0;
     virtual double GetDOutflowDp(double pressure, double z) const = 0;
     virtual double GetIntegralpD2S (double area, double z) const = 0;
     virtual void   SetPeriod(double period) = 0;
-    virtual void   SetAreas_and_length(double S_top, double S_bottom, double z) = 0;
+
+    // Technically speaking, this doesnt at all belong in "material"...
+    // but for now, we'll keep it here, since there was already 
+    // an odd coupling to the spatial characteristics.
+    void SetSpatialCharacteristics(cvOneD::SegmentSpatialCharacteristics const& in){ spatialCharacteristics = in; };
+
+    // Computes the initial radius at the z location based on the spatial
+    // characteristics. Again, doesn't really belong here...
+    double Getr1(double z) const{
+      // linearly interpolated r
+      return cvOneD::getInterpolatedRadius(z, spatialCharacteristics);
+    }
 
     double GetProfileExponent() const {return profile_exponent;}
     double GetDensity() const {return density;}
@@ -112,6 +125,8 @@ class cvOneDMaterial{
 
   protected:
 
+    cvOneD::SegmentSpatialCharacteristics spatialCharacteristics;
+
     double density;
     double dynamicViscosity;
     double kinematicViscosity;

Code/Source/cvOneDMaterialLinear.cxx

@@ -115,27 +115,12 @@ double cvOneDMaterialLinear::GetEHR(double z)const{
   return ehr;
 }
 
-void cvOneDMaterialLinear::SetAreas_and_length(double S_top,double S_bottom,double z){
-  Stop = S_top;//inlet
-  Sbot = S_bottom;//outlet
-  len = z;
-}
-
 double cvOneDMaterialLinear::GetS1(double z)const{
   double r    = Getr1(z);
   double area = r*r*M_PI;
   return area;
 }
 
-double cvOneDMaterialLinear::Getr1(double z)const{
-  // Linearly interpolated r
-  double r_top = sqrt(Stop/M_PI);
-  double r_bot = sqrt(Sbot/M_PI);
-  double r     = ((z-len)/(-len))*(r_top - r_bot) + r_bot;
-  
-  return r;
-}
-
 double cvOneDMaterialLinear::GetDS1Dz(double z)const{
   // Since vessel geometry will always be linear in
   // space, this is just a constant

Code/Source/cvOneDMaterialLinear.h

@@ -85,7 +85,6 @@ class cvOneDMaterialLinear:public cvOneDMaterial{
     double PP1_;
 
     double GetS1( double z) const;
-    double Getr1( double z) const;
     double GetDS1Dz( double z) const;
     double GetDr1Dz(double z) const;
 };

Code/Source/cvOneDMaterialOlufsen.cxx

@@ -151,16 +151,6 @@ double cvOneDMaterialOlufsen::GetS1(double z)const{
   return area;  // slightly increased pressure/decreased area
 }
 
-double cvOneDMaterialOlufsen::Getr1(double z)const{
-  // linearly interpolated r
-  double r_top=sqrt(Stop/PI);
-  double r_bot=sqrt(Sbot/PI);
-  double r=((z-len)/(-len))*(r_top - r_bot) + r_bot;
-
-
-  return r;
-}
-
 double cvOneDMaterialOlufsen::GetDS1Dz(double z)const{
   double drdz=GetDr1Dz(z) ;
   double dsdr= 2.0*PI*Getr1(z);