oomph::SolidHelpers Namespace Reference

Namespace for solid mechanics helper functions. More...

Functions
template<class ELEMENT >
void	doc_2D_principal_stress (DocInfo &doc_info, SolidMesh *mesh_pt)

Detailed Description

Namespace for solid mechanics helper functions.

/////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////

Function Documentation

doc_2D_principal_stress()

template<class ELEMENT >

void oomph::SolidHelpers::doc_2D_principal_stress	(	DocInfo &	doc_info,
		SolidMesh *	mesh_pt
	)

Document the principal stresses in a 2D SolidMesh pointed to by mesh_pt, in the directory specified by the DocInfo object, in a format that can be processed with tecplot macro.

    {
      // Output principal stress vectors at the centre of all elements
      std::ofstream pos_file;
      std::ofstream neg_file;
      std::ostringstream filename;
      filename << doc_info.directory() << "/pos_principal_stress"
               << doc_info.number() << ".dat";
      pos_file.open(filename.str().c_str());
      filename.str("");
      filename << doc_info.directory() << "/neg_principal_stress"
               << doc_info.number() << ".dat";
      neg_file.open(filename.str().c_str());
 
      // Write dummy data in both so there's at lest one zone in each
      pos_file << 0.0 << " " << 0.0 << " " << 0.0 << " " << 0.0 << " "
               << std::endl;
      neg_file << 0.0 << " " << 0.0 << " " << 0.0 << " " << 0.0 << " "
               << std::endl;
 
 
      Vector<double> s(2);
      Vector<double> x(2);
      s[0] = 0.0;
      s[1] = 0.0;
      unsigned n_solid_element = mesh_pt->nelement();
      for (unsigned e = 0; e < n_solid_element; e++)
      {
        ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt->element_pt(e));
 
        // Get principal stress
        DenseMatrix<double> principal_stress_vector(2);
        Vector<double> principal_stress(2);
        el_pt->get_principal_stress(
          s, principal_stress_vector, principal_stress);
 
        // Get position of centre of element
        el_pt->interpolated_x(s, x);
 
        // compute vectors at 45 degree for nearly hydrostatic pressure state
        DenseMatrix<double> rot(2);
 
        bool hydrostat = false;
 
        // Max. relative difference between principal stresses
        // required to classify stress state as non-hydrostatic: 1%
        double dev_max = 1.0e-2;
        if (principal_stress[0] != 0.0)
        {
          if (std::fabs((principal_stress[0] - principal_stress[1]) /
                        principal_stress[0]) < dev_max)
          {
            hydrostat = true;
            double Cos = cos(0.25 * 3.14159);
            double Sin = sin(0.25 * 3.14159);
            rot(0, 0) = Cos * principal_stress_vector(0, 0) -
                        Sin * principal_stress_vector(0, 1);
            rot(0, 1) = Sin * principal_stress_vector(0, 0) +
                        Cos * principal_stress_vector(0, 1);
            rot(1, 0) = Cos * principal_stress_vector(1, 0) -
                        Sin * principal_stress_vector(1, 1);
            rot(1, 1) = Sin * principal_stress_vector(1, 0) +
                        Cos * principal_stress_vector(1, 1);
          }
        }
 
        // Loop over two principal stresses:
        for (unsigned i = 0; i < 2; i++)
        {
          if (principal_stress[i] > 0.0)
          {
            pos_file << x[0] << " " << x[1] << " "
                     << principal_stress_vector(i, 0) << " "
                     << principal_stress_vector(i, 1) << std::endl;
            pos_file << x[0] << " " << x[1] << " "
                     << -principal_stress_vector(i, 0) << " "
                     << -principal_stress_vector(i, 1) << std::endl;
            if (hydrostat)
            {
              pos_file << x[0] << " " << x[1] << " " << rot(i, 0) << " "
                       << rot(i, 1) << std::endl;
              pos_file << x[0] << " " << x[1] << " " << -rot(i, 0) << " "
                       << -rot(i, 1) << std::endl;
            }
          }
          else
          {
            neg_file << x[0] << " " << x[1] << " "
                     << principal_stress_vector(i, 0) << " "
                     << principal_stress_vector(i, 1) << std::endl;
            neg_file << x[0] << " " << x[1] << " "
                     << -principal_stress_vector(i, 0) << " "
                     << -principal_stress_vector(i, 1) << std::endl;
            if (hydrostat)
            {
              neg_file << x[0] << " " << x[1] << " " << rot(i, 0) << " "
                       << rot(i, 1) << std::endl;
              neg_file << x[0] << " " << x[1] << " " << -rot(i, 0) << " "
                       << -rot(i, 1) << std::endl;
            }
          }
        }
      }
 
      pos_file.close();
      neg_file.close();
    }

References cos(), oomph::DocInfo::directory(), e(), oomph::Mesh::element_pt(), boost::multiprecision::fabs(), MergeRestartFiles::filename, i, oomph::Mesh::nelement(), oomph::DocInfo::number(), rot(), s, sin(), and plotDoE::x.