Fix matrix generation and shape factor #99

src/main.cpp

@@ -20,9 +20,9 @@ int main(int argc, char* argv[])
 
   bool WarmStartingFlag = false;
   int Resolution = 0;
-  double nu1 = 0.0, nu2 = 0.0, CompositeYoungs = 0.0, CompositePoissonsRatio = 0.0,
-         ShapeFactor = 0.0, ElasticComplianceCorrection = 0.0, GridSize = 0.0, Tolerance = 0.0,
-         E1 = 0.0, E2 = 0.0, LateralLength = 0.0, Delta = 0.0;
+  double nu1 = 0.0, nu2 = 0.0, CompositeYoungs = 0.0, ShapeFactor = 0.0,
+         ElasticComplianceCorrection = 0.0, GridSize = 0.0, Tolerance = 0.0, E1 = 0.0, E2 = 0.0,
+         LateralLength = 0.0, Delta = 0.0;
   bool RandomTopologyFlag = false;
   bool RandomSeedFlag = false;
   double Hurst = 0.0;
@@ -32,10 +32,10 @@ int main(int argc, char* argv[])
   int MaxIteration = 0;
   bool PressureGreenFunFlag = false;
 
-  MIRCO::SetParameters(E1, E2, LateralLength, nu1, nu2, CompositeYoungs, CompositePoissonsRatio,
-      ShapeFactor, ElasticComplianceCorrection, GridSize, Tolerance, Delta, TopologyFilePath,
-      Resolution, InitialTopologyStdDeviation, inputFileName, RandomTopologyFlag, Hurst,
-      RandomSeedFlag, RandomGeneratorSeed, WarmStartingFlag, MaxIteration, PressureGreenFunFlag);
+  MIRCO::SetParameters(E1, E2, LateralLength, nu1, nu2, CompositeYoungs, ShapeFactor,
+      ElasticComplianceCorrection, GridSize, Tolerance, Delta, TopologyFilePath, Resolution,
+      InitialTopologyStdDeviation, inputFileName, RandomTopologyFlag, Hurst, RandomSeedFlag,
+      RandomGeneratorSeed, WarmStartingFlag, MaxIteration, PressureGreenFunFlag);
 
   // Identical Vectors/Matricies, therefore only created one here.
   int ngrid = int(ceil((LateralLength - (GridSize / 2)) / GridSize));
@@ -60,8 +60,8 @@ int main(int argc, char* argv[])
   double pressure = 0.0;
 
   MIRCO::Evaluate(pressure, Delta, LateralLength, GridSize, Tolerance, MaxIteration,
-      CompositeYoungs, CompositePoissonsRatio, WarmStartingFlag, ElasticComplianceCorrection,
-      topology, max_and_mean.max_, meshgrid, PressureGreenFunFlag);
+      CompositeYoungs, WarmStartingFlag, ElasticComplianceCorrection, topology, max_and_mean.max_,
+      meshgrid, PressureGreenFunFlag);
 
   std::cout << "Mean pressure is: " << std::to_string(pressure) << std::endl;
 

src/mirco_evaluate.cpp

@@ -22,10 +22,9 @@
 
 
 void MIRCO::Evaluate(double& pressure, double Delta, double LateralLength, double GridSize,
-    double Tolerance, int MaxIteration, double CompositeYoungs, double CompositePoissonsRatio,
-    bool WarmStartingFlag, double ElasticComplianceCorrection,
-    Teuchos::SerialDenseMatrix<int, double>& topology, double zmax, std::vector<double>& meshgrid,
-    bool PressureGreenFunFlag)
+    double Tolerance, int MaxIteration, double CompositeYoungs, bool WarmStartingFlag,
+    double ElasticComplianceCorrection, Teuchos::SerialDenseMatrix<int, double>& topology,
+    double zmax, std::vector<double>& meshgrid, bool PressureGreenFunFlag)
 {
   // Initialise the area vector and force vector. Each element containing the
   // area and force calculated at every iteration.
@@ -71,8 +70,7 @@ void MIRCO::Evaluate(double& pressure, double Delta, double LateralLength, doubl
 
     // Construction of the Matrix A
     MIRCO::MatrixGeneration matrix1;
-    matrix1.SetUpMatrix(
-        A, xv0, yv0, GridSize, CompositeYoungs, CompositePoissonsRatio, n0, PressureGreenFunFlag);
+    matrix1.SetUpMatrix(A, xv0, yv0, GridSize, CompositeYoungs, n0, PressureGreenFunFlag);
 
     // Second predictor for contact set
     // @{

src/mirco_evaluate.h

@@ -16,7 +16,6 @@ namespace MIRCO
    * @param Tolerance Tolerance for the convergence of force
    * @param MaxIteration Maximum number of iterations for the force to converge
    * @param CompositeYoungs Composite Young's modulus
-   * @param CompositePoissonsRatio Composite Poisson's ratio
    * @param WarmStartingFlag Warm-Starter Flag
    * @param ElasticComplianceCorrection Elastic compliance correction
    * @param topology Topology matrix containing heights
@@ -26,10 +25,9 @@ namespace MIRCO
    * point force
    */
   void Evaluate(double& pressure, double Delta, double LateralLength, double GridSize,
-      double Tolerance, int MaxIteration, double CompositeYoungs, double CompositePoissonsRatio,
-      bool WarmStartingFlag, double ElasticComplianceCorrection,
-      Teuchos::SerialDenseMatrix<int, double>& topology, double zmax, std::vector<double>& meshgrid,
-      bool PressureGreenFunFlag);
+      double Tolerance, int MaxIteration, double CompositeYoungs, bool WarmStartingFlag,
+      double ElasticComplianceCorrection, Teuchos::SerialDenseMatrix<int, double>& topology,
+      double zmax, std::vector<double>& meshgrid, bool PressureGreenFunFlag);
 }  // namespace MIRCO
 
 

src/mirco_matrixsetup.cpp

@@ -7,7 +7,7 @@
 
 void MIRCO::MatrixGeneration::SetUpMatrix(Teuchos::SerialDenseMatrix<int, double>& A,
     std::vector<double> xv0, std::vector<double> yv0, double GridSize, double CompositeYoungs,
-    double CompositePoissonsRatio, int systemsize, bool PressureGreenFunFlag)
+    int systemsize, bool PressureGreenFunFlag)
 {
   double pi = M_PI;
   if (PressureGreenFunFlag)
@@ -27,7 +27,7 @@ void MIRCO::MatrixGeneration::SetUpMatrix(Teuchos::SerialDenseMatrix<int, double
                    n * log((sqrt(n * n + l * l) + l) / (sqrt(n * n + k * k) + k)));
       }
     }
-    A.scale((1 - pow(CompositePoissonsRatio, 2)) / (pi * CompositeYoungs));
+    A.scale(1 / (pi * CompositeYoungs));
   }
   else
   {

src/mirco_matrixsetup.h

@@ -18,14 +18,13 @@ namespace MIRCO
      * @param yv0 y-coordinates of the points in contact in the previous iteration.
      * @param GridSize Grid size (length of each cell)
      * @param CompositeYoungs The composite Young's modulus
-     * @param CompositePoissonsRatio The composite Poisson's ratio
      * @param systemsize Number of nodes predicted to be in contact
      * @param PressureGreenFunFlag Flag to use Green function based on uniform pressure instead of
      * point force
      */
     void SetUpMatrix(Teuchos::SerialDenseMatrix<int, double>& A, std::vector<double> xv0,
-        std::vector<double> yv0, double GridSize, double CompositeYoungs,
-        double CompositePoissonsRatio, int systemsize, bool PressureGreenFunFlag);
+        std::vector<double> yv0, double GridSize, double CompositeYoungs, int systemsize,
+        bool PressureGreenFunFlag);
     MatrixGeneration() = default;
   };
 }  // namespace MIRCO

src/mirco_setparameters.cpp

@@ -11,11 +11,11 @@
 #include "mirco_filesystem_utils.h"
 
 void MIRCO::SetParameters(double& E1, double& E2, double& LateralLength, double& nu1, double& nu2,
-    double& CompositeYoungs, double& CompositePoissonsRatio, double& ShapeFactor,
-    double& ElasticComplianceCorrection, double& GridSize, double& Tolerance, double& Delta,
-    std::string& TopologyFilePath, int& Resolution, double& InitialTopologyStdDeviation,
-    const std::string& inputFileName, bool& RandomTopologyFlag, double& Hurst, bool& RandomSeedFlag,
-    int& RandomGeneratorSeed, bool& WarmStartingFlag, int& MaxIteration, bool& PressureGreenFunFlag)
+    double& CompositeYoungs, double& ShapeFactor, double& ElasticComplianceCorrection,
+    double& GridSize, double& Tolerance, double& Delta, std::string& TopologyFilePath,
+    int& Resolution, double& InitialTopologyStdDeviation, const std::string& inputFileName,
+    bool& RandomTopologyFlag, double& Hurst, bool& RandomSeedFlag, int& RandomGeneratorSeed,
+    bool& WarmStartingFlag, int& MaxIteration, bool& PressureGreenFunFlag)
 {
   Teuchos::RCP<Teuchos::ParameterList> parameterList = Teuchos::rcp(new Teuchos::ParameterList());
   Teuchos::updateParametersFromXmlFile(inputFileName, parameterList.ptr());
@@ -44,14 +44,6 @@ void MIRCO::SetParameters(double& E1, double& E2, double& LateralLength, double&
   // Composite Young's modulus.
   CompositeYoungs = pow(((1 - pow(nu1, 2)) / E1 + (1 - pow(nu2, 2)) / E2), -1);
 
-  // Composite Shear modulus
-  double G1 = E1 / (2 * (1 + nu1));
-  double G2 = E2 / (2 * (1 + nu2));
-  double CompositeShear = pow(((2 - nu1) / (4 * G1) + (2 - nu2) / (4 * G2)), -1);
-
-  // Composite Poisson's ratio
-  CompositePoissonsRatio = CompositeYoungs / (2 * CompositeShear) - 1;
-
   // Shape factors (See section 3.3 of https://doi.org/10.1007/s00466-019-01791-3)
   // These are the shape factors to calculate the elastic compliance correction of the micro-scale
   // contact constitutive law for various resolutions.
@@ -63,7 +55,7 @@ void MIRCO::SetParameters(double& E1, double& E2, double& LateralLength, double&
   // http://dx.doi.org/10.1134/s1029959914040109
   const std::map<int, double> shape_factors_pressure{{1, 0.961389237917602}, {2, 0.924715342432435},
       {3, 0.899837531880697}, {4, 0.884976751041942}, {5, 0.876753783192863},
-      {6, 0.872397956576882}, {7, 0.871958228537090}, {8, 0.882669916668780}};
+      {6, 0.872397956576882}, {7, 0.871958228537090}, {8, 0.8689982669426167}};
 
   // The following force based constants are taken from Table 1 of Bonari et al. (2020).
   // https://doi.org/10.1007/s00466-019-01791-3

src/mirco_setparameters.h

@@ -16,7 +16,6 @@ namespace MIRCO
    * @param nu1 Poisson's ratio of body 1
    * @param nu2 Poisson's ratio of body 2
    * @param CompositeYoungs The composite Young's modulus
-   * @param CompositePoissonsRatio The composite Poisson's ratio
    * @param alpha Correction factor for this problem. Depends on the resolution.
    * @param ElasticComplianceCorrection Elastic compliance correction
    * @param GridSize Grid size (length of each cell)
@@ -40,12 +39,11 @@ namespace MIRCO
    * point force
    */
   void SetParameters(double& E1, double& E2, double& LateralLength, double& nu1, double& nu2,
-      double& CompositeYoungs, double& CompositePoissonsRatio, double& alpha,
-      double& ElasticComplianceCorrection, double& GridSize, double& Tolerance, double& Delta,
-      std::string& TopologyFilePath, int& Resolution, double& InitialTopologyStdDeviation,
-      const std::string& inputFileName, bool& RandomTopologyFlag, double& Hurst,
-      bool& RandomSeedFlag, int& RandomGeneratorSeed, bool& WarmStartingFlag, int& MaxIteration,
-      bool& PressureGreenFunFlag);
+      double& CompositeYoungs, double& alpha, double& ElasticComplianceCorrection, double& GridSize,
+      double& Tolerance, double& Delta, std::string& TopologyFilePath, int& Resolution,
+      double& InitialTopologyStdDeviation, const std::string& inputFileName,
+      bool& RandomTopologyFlag, double& Hurst, bool& RandomSeedFlag, int& RandomGeneratorSeed,
+      bool& WarmStartingFlag, int& MaxIteration, bool& PressureGreenFunFlag);
 }  // namespace MIRCO
 
 #endif  // SRC_SETPARAMETERS_H_