matlabElasticityGraphics.hh

Go to the documentation of this file.

 
#ifndef MATLAB_ELASTICITY_HH_
#define MATLAB_ELASTICITY_HH_
 
#include "graphics.hh"
#include <iostream>
 
using concepts::Real;
 
namespace graphics {
 
class MatlabElasticityGraphics: public MatlabBinaryGraphics {
 
public:
 
  template<class G, class F>
  MatlabElasticityGraphics(const concepts::Space<G>& spc,
      const std::string& filename, const concepts::Vector<F>& sol_u1,
      const concepts::Vector<F>& sol_u2, F lambda, F mu,
      bool plain_stress);
 
  virtual ~MatlabElasticityGraphics() {
  }
  ;
 
};
 
template<class G, class F>
MatlabElasticityGraphics::MatlabElasticityGraphics(
    const concepts::Space<G>& spc,
    const std::string& filename,
    const concepts::Vector<F>& sol_u1,
    const concepts::Vector<F>& sol_u2,
    F lambda, F mu,
    bool plain_stress) :
    MatlabBinaryGraphics(filename) {
 
  //read necessary informations
  hp2D::Grad<F> grad;
  DenseMatrixCollection<G> dmc = this->getSpace(spc);
  //std::map<std::string, concepts::ElementMatrix<G>*> map = dmc.getMapping();
 
  concepts::RCP<concepts::ElementMatrix<F> > u1 = this->getSolution(spc, sol_u1);
  concepts::RCP<concepts::ElementMatrix<F> > u2 = this->getSolution(spc, sol_u2);
  concepts::RCP<concepts::ElementMatrix<F> > grad_u1 = this->getSolution(spc, sol_u1, &grad);
  concepts::RCP<concepts::ElementMatrix<F> > grad_u2 = this->getSolution(spc, sol_u2, &grad);
 
//  concepts::ElementMatrix<F> u1 = *getSolution(spc, sol_u1);
//  concepts::ElementMatrix<F> u2 = *getSolution(spc, sol_u2);
//
//  concepts::ElementMatrix<F> grad_u1 = this->getSolution(spc, sol_u1, &grad);
//  concepts::ElementMatrix<F> grad_u2 = this->getSolution(spc, sol_u2, &grad);
 
  uint n = grad_u1->n();
 
  // compute strains and stresses
  concepts::Vector<F> strain11_(n);
  concepts::Vector<F> strain22_(n);
  concepts::Vector<F> strain12_(n);
  concepts::Vector<F> stress11_(n);
  concepts::Vector<F> stress22_(n);
  concepts::Vector<F> stress12_(n);
  concepts::Vector<F> stress33_(n);
 
  for (uint i = 0; i < n; i++) {
    strain11_(i) = grad_u1->operator ()(0, i);
    strain22_(i) = grad_u2->operator ()(1, i);
    strain12_(i) = 0.5
        * (grad_u2->operator ()(0, i) + grad_u1->operator ()(1, i));
  }
  if (plain_stress)
    for(uint i = 0; i < n; i++) {
      stress11_(i) = 2*mu*(lambda/(lambda+2*mu) + 1) * strain11_(i)
                          + 2*mu*lambda/(lambda+2*mu) * strain22_(i);
      stress22_(i) = 2*mu*lambda/(lambda+2*mu) * strain11_(i)
                          + 2*mu*(lambda/(lambda+2*mu) + 1) * strain22_(i);
      stress12_(i) = 2 * mu * strain12_(i);
    }
  else
    for (uint i = 0; i < n; i++) {
      stress11_(i) = (lambda + 2 * mu) * strain11_(i) + lambda * strain22_(i);
      stress22_(i) = lambda * strain11_(i) + (lambda + 2 * mu) * strain22_(i);
      stress12_(i) = 2 * mu * strain12_(i);
      stress33_(i) = lambda * (strain11_(i) + strain22_(i) );
    }
 
 
  // convert u1, u2 to vectors
  concepts::Vector<F> u1_(n), u2_(n);
  for (uint i = 0; i < n; i++){
    u1_(i) = u1->operator()(0,i);
    u2_(i) = u2->operator()(0,i);
  }
 
 
  // compute the deformed mesh
//  concepts::ElementMatrix<F> x_new = *map["x"];
//  x_new.add(u1);
//  concepts::ElementMatrix<F> y_new = *map["y"];
//  y_new.add(u1);
 
 
 
  // write data to output
//  this->storeData_(true, dmc->getMapping());
  this->storeData_(dmc);
 
//  this->add(x_new, "x_new");
//  this->add(y_new, "y_new");
  this->add(stress11_, "stress11");
  this->add(stress22_, "stress22");
  this->add(stress12_, "stress12");
  if (!plain_stress){
    this->add(stress33_, "stress33");
  }
  this->add(u1_, "displacement_u1");
  this->add(u2_, "displacement_u2");
 
}
 
} // namespace graphics
 
#endif