normmat.cpp

Finite element program for mechanical problem. It can solve various problem in solid problem

                                  6*k*epschr1*


                                    pow(epsc1, 


                                      2)*(as2*ch*e0c*epschr1 + 


                                        d*(n0 - 


                                        e0c*(as1*epsn + 


                                        as2*(epschr1 + epsn))))))))) - 


      6*as2*e0c*(-2*ch + h)*(epschr1 - (ch*epschr1)/d + 


            epsn + (ch*(6*as2*ch*e0c*k*epschr1*pow(epsc1, 2) + 


                      3*d*k*pow(epsc1, 


                          2)*(n0 + as1*e0c*(epschr1 - epsn) - 


                            as2*e0c*(epschr1 + epsn)) + 


                      sqrt(3.)*


                        sqrt(k*


                            pow(epsc1, 


                              2)*(3*k*pow(epsc1, 2)*


                                  pow(-2*as2*ch*e0c*epschr1 + 


                                      d*(-n0 + 


                                        e0c*(as1*(-epschr1 + epsn) + 


                                        as2*(epschr1 + epsn))), 2) - 


                                2*(as2*ch + as1*d)*


                                  e0c*(b*pow(d, 2)*fc*


                                        pow(epschr1 + epsn, 


                                        2)*(4*(-1 + k)*epschr1 + (6 - 9*k)*


                                        epsc1 + 4*(-1 + k)*epsn) + 


                                      6*k*epschr1*


                                        pow(epsc1, 


                                        2)*(as2*ch*e0c*epschr1 + 


                                        d*(n0 - e0c*(as1*epsn + 


                                        as2*(epschr1 + epsn)))))))))/(6.*


                  d*(as2*ch + as1*d)*e0c*k*pow(epsc1, 2))) - (b*d*fc*


            pow(epschr1 + epsn, 


              2)*(2*d*(epschr1 + 


                      epsn)*((-1 + k)*epschr1 + (2 - 3*k)*epsc1 + (-1 + k)*


                        epsn) - 


                h*(4*(-1 + k)*epschr1 + (6 - 9*k)*epsc1 + 


                      4*(-1 + k)*


                        epsn)*(epschr1 - (6*as2*ch*e0c*k*epschr1*


                              pow(epsc1, 2) + 


                            3*d*k*pow(epsc1, 


                                2)*(n0 + as1*e0c*(epschr1 - epsn) - 


                                  as2*e0c*(epschr1 + epsn)) + 


                            sqrt(3.)*


                              sqrt(k*


                                  pow(epsc1, 


                                    2)*(3*k*pow(epsc1, 2)*


                                        pow(-2*as2*ch*e0c*epschr1 + 


                                        d*(-n0 + e0c*(as1*(-epschr1 + epsn) + 


                                        as2*(epschr1 + epsn))), 2) - 


                                      2*(as2*ch + as1*d)*


                                        e0c*(b*pow(d, 2)*fc*


                                        pow(epschr1 + epsn, 


                                        2)*(4*(-1 + k)*epschr1 + (6 - 9*k)*


                                        epsc1 + 4*(-1 + k)*epsn) + 


                                        6*k*epschr1*


                                        pow(epsc1, 


                                        2)*(as2*ch*e0c*epschr1 + 


                                        d*(n0 - e0c*(as1*epsn + 


                                        as2*(epschr1 + epsn))))))))/(6.*(as2*


                                  ch + as1*d)*e0c*k*pow(epsc1, 2)))))/(k*


            pow(epsc1, 2)*


            pow(epschr1 - 


                1/(6.*(as2*ch + as1*d)*e0c*k*pow(epsc1, 2))*(6*as2*ch*e0c*k*


                        epschr1*pow(epsc1, 2) + 


                      3*d*k*pow(epsc1, 


                          2)*(n0 + as1*e0c*(epschr1 - epsn) - 


                            as2*e0c*(epschr1 + epsn)) + 


                      sqrt(3.)*


                        sqrt(k*


                            pow(epsc1, 


                              2)*(3*k*pow(epsc1, 2)*


                                  pow(-2*as2*ch*e0c*epschr1 + 


                                      d*(-n0 + 


                                        e0c*(as1*(-epschr1 + epsn) + 


                                        as2*(epschr1 + epsn))), 2) - 


                                2*(as2*ch + as1*d)*


                                  e0c*(b*pow(d, 2)*fc*


                                        pow(epschr1 + epsn, 


                                        2)*(4*(-1 + k)*epschr1 + (6 - 9*k)*


                                        epsc1 + 4*(-1 + k)*epsn) + 


                                      6*k*epschr1*


                                        pow(epsc1, 


                                        2)*(as2*ch*e0c*epschr1 + 


                                        d*(n0 - e0c*(as1*epsn + 


                                        as2*(epschr1 + epsn)))))))), 2))) + e0c*(-(as1*(-cd + h/2.)) + as2*(-ch + h/2.))*epsn;





  }


  else{


    mhr1=1000000000.0;


	 }


  return mhr1;


  }


  


 //TLACENA VYZTUZ + NORMALOVA SILA + TECENI V As1 + x < h


  


 //pretvoreni betonu, kdy dochazi k teceni vyztuze As1, ale netece As2 a ani nedochazi k drceni betonu


double normmat::compute_epschr2 ()


{


  double disk1,epsct,disk2,epschr2;





  disk1 =e0c*(3*e0c*


        pow(k*pow(epsc1, 


                2)*(-2*as2*ch*fy - as1*d*fy + as2*d*fy + d*n0 + 


                  2*as2*(ch - d)*e0c*epsn) + 


            b*pow(d, 2)*


              fc*(-2*(-1 + k)*pow(epsc1, 2) + (2 - 3*k)*epsc1*epsn + 


                  2*(-1 + k)*pow(epsn, 2)), 


          2) - (2*as2*(-ch + d)*e0c*k*pow(epsc1, 2) - 


            b*pow(d, 2)*fc*((-2 + k)*epsc1 + 4*(-1 + k)*epsn))*(-6*as2*k*


              pow(epsc1, 2)*(fy - e0c*epsn)*(ch*fy + (-ch + d)*e0c*epsn) + 


            d*(6*fy*k*n0*pow(epsc1, 2) - 


                  6*as1*fy*k*pow(epsc1, 2)*(fy - e0c*epsn) + 


                  e0c*(-6*k*n0*pow(epsc1, 2)*epsn + 


                        b*d*fc*(-4*(-1 + k)*pow(epsc1, 3) + 


                              3*(-2 + 3*k)*epsc1*pow(epsn, 2) - 


                              4*(-1 + k)*pow(epsn, 3))))));


  


  if (mhr1 == 1000000000.0 || disk1 < 0.0){





    epsct=-1.0;





  }





  else {





    epsct=(-3*e0c*k*pow(epsc1,2) * (-2*as2*ch*fy - as1*d*fy + as2*d*fy + d*n0 + 2*as2*(ch - d)*e0c*epsn) + 


	   3*b*pow(d, 2)*e0c*fc*(2*(-1 + k)*pow(epsc1, 2) + (-2 + 3*k)*epsc1*epsn - 2*(-1 + k)*pow(epsn, 2)) + 


	   sqrt(3.)*sqrt(e0c*(3*e0c*pow(k*pow(epsc1,2)*(-2*as2*ch*fy - as1*d*fy + as2*d*fy + d*n0 + 


                              2*as2*(ch - d)*e0c*epsn) + 


                        b*pow(d, 2)*


                          fc*(-2*(-1 + k)*pow(epsc1, 2) + (2 - 3*k)*epsc1*


                                epsn + 2*(-1 + k)*pow(epsn, 2)), 


                      2) - (2*as2*(-ch + d)*e0c*k*pow(epsc1, 2) - 


                        b*pow(d, 2)*


                          fc*((-2 + k)*epsc1 + 4*(-1 + k)*epsn))*(-6*as2*k*


                          pow(epsc1, 


                            2)*(fy - e0c*epsn)*(ch*fy + (-ch + d)*e0c*epsn) + 


                        d*(6*fy*k*n0*pow(epsc1, 2) - 


                              6*as1*fy*k*pow(epsc1, 2)*(fy - e0c*epsn) + 


                              e0c*(-6*k*n0*pow(epsc1, 2)*epsn + 


                                    b*d*fc*(-4*(-1 + k)*pow(epsc1, 3) + 


                                        3*(-2 + 3*k)*epsc1*pow(epsn, 2) - 


					    4*(-1 + k)*pow(epsn, 3))))))))/


      (3.*e0c*(2*as2*(ch - d)*e0c*k*pow(epsc1, 2) + b*pow(d, 2)*fc*((-2 + k)*epsc1 + 4*(-1 + k)*epsn)));





  }





  disk2 =e0c*(3*e0c*


        pow(as2*k*


              pow(epsc1, 2)*(-2*ch*fy + d*fy + 2*ch*e0c*epsn - 2*d*e0c*epsn) +


             d*(-(as1*fy*k*pow(epsc1, 2)) + 


                  2*b*d*fc*(pow(epsct, 2) + (epsc1 - epsn)*epsn) + 


                  k*(n0*pow(epsc1, 2) + 


                        b*d*fc*(-2*pow(epsct, 2) - 3*epsc1*epsn + 


                              2*pow(epsn, 2)))), 


          2) - (2*as2*(-ch + d)*e0c*k*pow(epsc1, 2) + 


            b*pow(d, 2)*


              fc*((-2 + 3*k)*epsc1 - 4*(-1 + k)*(epsct + epsn)))*(-6*as2*k*


              pow(epsc1, 2)*(fy - e0c*epsn)*(ch*fy + (-ch + d)*e0c*epsn) + 


            d*(6*fy*k*n0*pow(epsc1, 2) - 


                  6*as1*fy*k*pow(epsc1, 2)*(fy - e0c*epsn) + 


                  e0c*(-6*k*n0*pow(epsc1, 2)*epsn + 


                        b*d*fc*(-4*(-1 + k)*pow(epsct, 3) + 


                              3*(-2 + 3*k)*epsc1*pow(epsn, 2) - 


                              4*(-1 + k)*pow(epsn, 3))))));





  if (mhr1 == 1000000000.0 || disk2 < 0.0 || epsct == -1.0){





	  epschr2=-1.0;


  }


  else{





    epschr2=(-3*as2*e0c*k*


          pow(epsc1, 2)*(-2*ch*fy + d*fy + 2*ch*e0c*epsn - 2*d*e0c*epsn) + 


        3*d*e0c*(k*(as1*fy - n0)*pow(epsc1, 2) - 2*b*d*fc*pow(epsct, 2) + 


              2*b*d*fc*k*pow(epsct, 2) + b*d*fc*(-2 + 3*k)*epsc1*epsn - 


              2*b*d*fc*(-1 + k)*pow(epsn, 2)) + 


        sqrt(3.)*sqrt(e0c*(3*e0c*


                    pow(as2*k*


                          pow(epsc1, 


                            2)*(-2*ch*fy + d*fy + 2*ch*e0c*epsn - 


                              2*d*e0c*epsn) + 


                        d*(-as1*fy*k*pow(epsc1, 2) + 


                              2*b*d*fc*(pow(epsct, 2) + (epsc1 - epsn)*epsn) +


                               k*(n0*pow(epsc1, 2) + 


                                    b*d*fc*(-2*pow(epsct, 2) - 3*epsc1*epsn + 


                                        2*pow(epsn, 2)))), 


                      2) - (2*as2*(-ch + d)*e0c*k*pow(epsc1, 2) + 


                        b*pow(d, 2)*


                          fc*((-2 + 3*k)*epsc1 - 


                              4*(-1 + k)*(epsct + epsn)))*(-6*as2*k*


                          pow(epsc1, 


                            2)*(fy - e0c*epsn)*(ch*fy + (-ch + d)*e0c*epsn) + 


                        d*(6*fy*k*n0*pow(epsc1, 2) - 


                              6*as1*fy*k*pow(epsc1, 2)*(fy - e0c*epsn) + 


                              e0c*(-6*k*n0*pow(epsc1, 2)*epsn + 


                                    b*d*fc*(-4*(-1 + k)*pow(epsct, 3) + 


                                        3*(-2 + 3*k)*epsc1*pow(epsn, 2) - 


                                        4*(-1 + k)*pow(epsn, 3))))))))/(3.*


        e0c*(2*as2*(ch - d)*e0c*k*pow(epsc1, 2) + 


	     b*pow(d, 2)*fc*((2 - 3*k)*epsc1 + 4*(-1 + k)*(epsct + epsn))));


  }


  return epschr2;


}


 


 


 //hodnota ohyboveho momentu, kdy dochazi k teceni vyztuze As1, ale netece As2 a ani nedochazi k drceni betonu 


double normmat::compute_mhr2 ()


{


	double mhr2;	


	


	if (epschr2 < 0.0 && epschr2!=-1.0){





    mhr2=0.08333333333333333*(6*as1*


        e0c*(-2*cd + h)*(epsn + 


            1/(6.*(as2*ch + as1*d)*e0c*k*pow(epsc1, 2))*(6*as2*ch*e0c*k*


                    epschr2*pow(epsc1, 2) + 


                  3*d*k*pow(epsc1, 


                      2)*(n0 + as1*e0c*(epschr2 - epsn) - 


                        as2*e0c*(epschr2 + epsn)) + 


                  sqrt(3.)*


                    sqrt(k*pow(epsc1, 


                          2)*(3*k*pow(epsc1, 2)*


                              pow(-2*as2*ch*e0c*epschr2 + 


                                  d*(-n0 + 


                                        e0c*(as1*(-epschr2 + epsn) + 


                                        as2*(epschr2 + epsn))), 2) - 


                            2*(as2*ch + as1*d)*


                              e0c*(b*pow(d, 2)*fc*


                                    pow(epschr2 + epsn, 


                                      2)*(4*(-1 + k)*epschr2 + (6 - 9*k)*


                                        epsc1 + 4*(-1 + k)*epsn) + 


                                  6*k*epschr2*


                                    pow(epsc1, 


                                      2)*(as2*ch*e0c*epschr2 + 


                                        d*(n0 - 


                                        e0c*(as1*epsn + 


                                        as2*(epschr2 + epsn))))))))) - 


      6*as2*e0c*(-2*ch + h)*(epschr2 - (ch*epschr2)/d + 


            epsn + (ch*(6*as2*ch*e0c*k*epschr2*pow(epsc1, 2) + 


                      3*d*k*pow(epsc1, 


                          2)*(n0 + as1*e0c*(epschr2 - epsn) - 


                            as2*e0c*(epschr2 + epsn)) + 


                      sqrt(3.)*


                        sqrt(k*


                            pow(epsc1, 


                              2)*(3*k*pow(epsc1, 2)*


                                  pow(-2*as2*ch*e0c*epschr2 + 


                                      d*(-n0 + 


                                        e0c*(as1*(-epschr2 + epsn) + 


                                        as2*(epschr2 + epsn))), 2) - 


                                2*(as2*ch + as1*d)*


                                  e0c*(b*pow(d, 2)*fc*


                                        pow(epschr2 + epsn, 


                                        2)*(4*(-1 + k)*epschr2 + (6 - 9*k)*


                                        epsc1 + 4*(-1 + k)*epsn) + 


                                      6*k*epschr2*


                                        pow(epsc1, 


                                        2)*(as2*ch*e0c*epschr2 + 


                                        d*(n0 - e0c*(as1*epsn + 


                                        as2*(epschr2 + epsn)))))))))/(6.*


                  d*(as2*ch + as1*d)*e0c*k*pow(epsc1, 2))) - (b*d*fc*


            pow(epschr2 + epsn, 


              2)*(2*d*(epschr2 + 


                      epsn)*((-1 + k)*epschr2 + (2 - 3*k)*epsc1 + (-1 + k)*


                        epsn) - 


                h*(4*(-1 + k)*epschr2 + (6 - 9*k)*epsc1 + 


                      4*(-1 + k)*


                        epsn)*(epschr2 - (6*as2*ch*e0c*k*epschr2*


                              pow(epsc1, 2) + 


                            3*d*k*pow(epsc1, 


                                2)*(n0 + as1*e0c*(epschr2 - epsn) - 


                                  as2*e0c*(epschr2 + epsn)) + 


                            sqrt(3.)*