[Getfem-users] 2d plane strain/stress analysis

[Umut Tabak]

Dear all,

I was wondering when using the 
asm_stiffness_matrix_for_linear_elasticity function for a 2d problem(on 
a mesh that is imported from gmsh), I would like to define a 2d 
vectorial problem(a plane strain or stress one actually). So I set the 
dimension to 2 since this is a 2s mesh. But when using this method in 
the format given in the format of the manual

invalid data mesh fem (Qdim=1 required)

It has been sometime I have not been using the library, this can be a 
simple mistake. The second question is that the thickness of the domain 
in say z direction, is that taken as unit when computing the stiffness 
matrices? And if I would like to define a different thickness, should I 
use generic assembly procedure for the beam bending problems. I guess 
there I can supply this as a parameter to the integration process.

Simple code is attached.

Best regards,
Umut


//
// standard headers
//
#include <iostream>
#include <vector>
#include <string>
#include <cstdlib>
#include <cassert>
#include <stdexcept>
//
// gmm++ headers, pay attention that this is a template
// library, so that no build is required.
//
#include <gmm/gmm.h>
#include <gmm/gmm_inoutput.h>
//
// getFem++ headers 
//
#include <getfem/getfem_import.h>
#include <getfem/dal_bit_vector.h>
#include <getfem/getfem_mesh.h>
#include <getfem/getfem_mesh_fem.h>
#include <getfem/getfem_mesh_im.h>
#include <getfem/getfem_integration.h>
#include <getfem/getfem_assembling.h>
#include <getfem/getfem_modeling.h>
#include <getfem/getfem_regular_meshes.h>
typedef getfem::modeling_standard_sparse_matrix sparse_matrix;
using bgeot::scalar_type;
using bgeot::size_type;
//using namespace getfem;
//using namespace gmm;
using namespace std;
//
int main(int argc, char** argv)
{ 
  try{
    getfem::mesh mshGmsh;
    import_mesh("gmshv2:./beamModel2d.msh", mshGmsh);
    getfem::mesh_fem mf1(mshGmsh);
    size_type dim = mshGmsh.dim();
    std::cout << "Mesh dimension:" << dim << endl;
    mf1.set_qdim(dim);
    mf1.set_finite_element(mshGmsh.convex_index(), 
                                     getfem::fem_descriptor("FEM_QK(2, 1)"));

    // check mesh dim and vectorial dimension equality
    if(mshGmsh.dim() == mf1.get_qdim())
    {
      cout << "dimensions are equal" << endl;
      cout << "number of dofs of the problem" << mf1.nb_dof() << endl;
    }
    // set integration methods
    getfem::mesh_im intMethods(mshGmsh);
    getfem::pintegration_method ppi = 
getfem::int_method_descriptor("IM_QUAD(3)");
    string im_name = getfem::name_of_int_method(ppi);
    cout << im_name << endl;
    intMethods.set_integration_method(ppi);
    // define the matrix to use in the assembly process of
    sparse_matrix SM(mf1.nb_dof(), mf1.nb_dof());
    sparse_matrix MM(mf1.nb_dof(), mf1.nb_dof());
    // define elastic constants
    vector<double> lambda(mf1.nb_dof(), 70e9*0.33/((1+0.33)*(1-2*0.33)));
    vector<double> mu(mf1.nb_dof(), 70e9/(2*(1+0.33)));
    // assembly process
    //asm_stiffness_matrix_for_linear_elasticity(SM, intMethods, mf1, mf1, 
lambda, mu);
    std::cout<<" ---  stiffness matrix assembling --- " << std::endl;
    //assem.assembly();
    getfem::asm_stiffness_matrix_for_linear_elasticity(SM, intMethods, mf1, 
mf1, lambda, mu);
    std::cout<<" ---  stiffness matrix assembled  --- " << std::endl;
    // std::cout<<" ---  stiffness matrix:  " << endl << SM << std::endl;
    //std::cout<<" ---       mass matrix assembling --- " << std::endl;
    getfem::asm_mass_matrix(MM, intMethods, mf1, mf1);
    std::cout<<" ---       mass matrix assembled  --- " << std::endl;
    // write out the matrices in Harwell-Boeing format
    string kFileName("Kmat.mm");
    string mFileName("Mmat.mm");
    // Harwell_Boeing_save(mFileName.c_str(), MM);  
    // Harwell_Boeing_save(kFileName.c_str(), SM);  
    // write out the matrices in Matrix-Market format
    gmm::csc_matrix<double> SM2,MM2;
    //gmm::copy(SM, SM2);
    gmm::copy(MM, MM2);
    MatrixMarket_save(mFileName.c_str(), MM2);
    //MatrixMarket_save(kFileName.c_str(), SM2);
  }
  catch(std::exception &e){
    std::cout << e.what() << std::endl;
  }
  return EXIT_SUCCESS;
}