Update Compute_MuGamma for parametrized_Laminate

inputs/computeMuGamma.parset

@@ -29,8 +29,7 @@
 
 nElements = 32 32
 
-gamma=100
-
+gamma=1000000
 
 #############################################
 #  Material parameters
@@ -49,6 +48,7 @@ lambda1=0.0
 ####  material_implementation("analytical_Example") ? 
 
 material_prestrain_imp= "parametrized_Laminate"
+
 #material_prestrain_imp= "analytical_Example"
 #material_prestrain_imp ="isotropic_bilayer"
 
@@ -60,13 +60,13 @@ material_prestrain_imp= "parametrized_Laminate"
 #  Prestrain parameters
 #############################################
 
-write_prestrainFunctions = true  # VTK norm of B ,
+write_prestrainFunctions = false  # VTK norm of B ,
 
 
 rho1=1.0
 alpha=2.0
 
-theta=0.541859296482412
+theta=0.125
 
 #theta = 0.3    # volume fraction   #default = 1.0/4.0
 #theta = 0.25   # volume fraction

src/Compute_MuGamma.cc

@@ -748,7 +748,7 @@ int main(int argc, char *argv[])
     log << "material_prestrain used: "<< imp  << std::endl;
     double gamma   = parameterSet.get<double>("gamma", 1.0);   
     double theta   = parameterSet.get<double>("theta", 1.0/4.0); 
-    double mu1     = parameterSet.get<double>("mu1", 10.0);
+    double mu1     = parameterSet.get<double>("mu1", 1.0);
     double beta    = parameterSet.get<double>("beta", 2.0); 
     double mu2     = beta*mu1;
     std::cout << "Gamma:" << gamma << std::endl;
@@ -887,25 +887,49 @@ int main(int argc, char *argv[])
 //     Vector zeroVec(basis.size());
 //     zeroVec = 0;
 //     std::cout << "TEST computeMuGamma:" << computeMuGamma(basis, zeroVec, gamma, muLocal, ForceLocal)<< std::endl;
-
-    std::cout << " --- print analytic solutions --- " << std::endl;
+    std::cout << " --- print analytic solutions(if possible) --- " << std::endl;
+    if (imp == "analytical_Example")
+    {
+        std::cout<< "analytical_Example" << std::endl;
+        double q1 = ((mu1*mu2)/6.0)/(theta*mu1+ (1.0- theta)*mu2);
+        double q2 = ((1.0-theta)*mu1+theta*mu2)/6.0;
+        double hm = mu1*(beta/(theta+(1-theta)*beta)) *(1.0/6.0);  
+        double am = mu1*((1-theta)+theta*beta)     *(1.0/6.0);
+        std::cout << "q1 : "     << q1 << std::endl;
+        std::cout << "q2 : "     << q2 << std::endl;
+        std::cout << "q3 should be between q1 and q2"  << std::endl;
+        std::cout << "hm : "     << hm << std::endl;
+        std::cout << "am : "     << am << std::endl;
+        if(mu_gamma > q2)
+        {
+            std::cout << "mu_gamma > q2!!.. (39) not satisfied" << std::endl;
+        }
+    }
+    if (imp == "parametrized_Laminate")
+    {
+        std::cout<< "parametrized_Laminate" << std::endl;
+        double hm = mu1*(beta/(theta+(1-theta)*beta));  
+        double am = mu1*((1-theta)+theta*beta);
+        double q1 = (1.0/6.0)*hm;
+        double q2 = (1.0/6.0)*am;
+        std::cout << "q1 : "     << q1 << std::endl;
+        std::cout << "q2 : "     << q2 << std::endl;
+        std::cout << "q3 should be between q1 and q2"  << std::endl;
+        std::cout << "hm : "     << hm << std::endl;
+        std::cout << "am : "     << am << std::endl;
+        if(mu_gamma > q2)
+        {
+            std::cout << "mu_gamma > q2!!.. (39) not satisfied" << std::endl;
+        }
+    }
     
-    double q1 = ((mu1*mu2)/6.0)/(theta*mu1+ (1.0- theta)*mu2);
-    double q2 = ((1.0-theta)*mu1+theta*mu2)/6.0;
-    double hm = mu1*(beta/(theta+(1-theta)*beta)) *(1.0/6.0);  
-    double am = mu1*((1-theta)+theta*beta)     *(1.0/6.0);
     
-    std::cout << "q1 : "     << q1 << std::endl;
-    std::cout << "q2 : "     << q2 << std::endl;
-    std::cout << "q3 should be between q1 and q2"  << std::endl;
-    std::cout << "hm : "     << hm << std::endl;
-    std::cout << "am : "     << am << std::endl;
     
+
+
     
-    if(mu_gamma > q2)
-    {
-        std::cout << "mu_gamma > q2!!.. (39) not satisfied" << std::endl;
-    }
+    
+
 
     
     std::cout << "beta : "     << beta << std::endl;