[Getfem-users] plate elements

[Umut Tabak]

Dear a users,
I tried to use the plate elements in getfem. I have simple code which 
seems to work I can get the dof counts right, as a dynamics person, I 
need the system matrices K and M, so can someone put some light on the 
extraction of system matrices for the below code snippet. Maybe the code 
could be buggy, but I should get the matrices from out of the model 
bricks I guess. Any help is appreciated.

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_partial_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>
#include <getfem/getfem_linearized_plates.h>
//
typedef getfem::modeling_standard_sparse_matrix sparse_matrix;
typedef getfem::modeling_standard_plain_vector plain_vector;
//using namespace getfem;
//using namespace gmm;
using namespace std;
using bgeot::scalar_type;
using bgeot::size_type;   /* = unsigned long */
//
int main(int argc, char** argv) { 
 // try {
   // if(argc!=3){
   //   throw std::runtime_error("Error in input parameters\nUsage: <binary> 
<filename.msh> <gmshFormatSpecifier>");
   // }
   // mesh definition
   getfem::mesh mesh;
   // integration specifications
        getfem::mesh_im mim(mesh), mim_subint(mesh);
        // mesh fem definitions for the fields on the plate
        // u3 out of plane definition
   // ut in plane, membrane definitions
   // theta rotation definitions
        getfem::mesh_fem mf_u3(mesh);
        getfem::mesh_fem mf_ut(mesh);
        getfem::mesh_fem mf_theta(mesh);
   // mesh fem for rhs and lame constants
        getfem::mesh_fem mf_rhs(mesh);
        getfem::mesh_fem mf_lame(mesh);
        // plate thickness
        scalar_type h=0.005;
   scalar_type lambda = 70e9*.33/((1+0.33)*(1-2*0.33));
   scalar_type mu     = 70e9/(2*(1+0.33));
   // mitc
        bool mitc = true;
        // eta for Kirchoff-Love '0' or Mindlin 'small'
        scalar_type eta = 0;

   //import_mesh(argv[1], argv[2], mshGmsh);
   // there is one mesh definition and different mesh_fem definitions
   // for different regions in definition
   // getfem::import_mesh("gmshv2:./platecoarse.msh", mesh);
   getfem::import_mesh("gmshv2:./shellMesh.msh", mesh);
   if( mesh.has_region(1)){
     std::cout << "Mesh has the specified region definitions" << std::endl;
   }
   // optional print out for the convexes(elements) of a specific
   // region
   cout << "Number of elements : " << mesh.nb_convex() << endl;
   // assign fem definitions to the regions
   mf_ut.set_qdim(2);
        mf_theta.set_qdim(2);
   // set the finite element method
        getfem::pfem pf_ut = getfem::fem_descriptor("FEM_QK(2,1)");
        getfem::pfem pf_u3 = getfem::fem_descriptor("FEM_QK(2,1)");
   getfem::pfem pf_theta = getfem::fem_descriptor("FEM_QK(2,1)");
        
        getfem::pintegration_method ppi = 
getfem::int_method_descriptor("IM_GAUSS_PARALLELEPIPED(2, 10)");
        getfem::pintegration_method ppi_ct = 
getfem::int_method_descriptor("IM_GAUSS_PARALLELEPIPED(2, 1)");

        mim.set_integration_method(mesh.convex_index(), ppi);
        mim_subint.set_integration_method(mesh.convex_index(), ppi_ct);
        mf_ut.set_finite_element(mesh.convex_index(), pf_ut);
        mf_u3.set_finite_element(mesh.convex_index(), pf_u3);
   mf_theta.set_finite_element(mesh.convex_index(), pf_theta);

   //
        mf_rhs.set_finite_element(mesh.convex_index(),pf_ut);

   //for(getfem::mr_visitor i(mshGmsh.region(1)); !i.finished(); ++i){
   //  lameConsts.set_finite_element(i.cv(),
        //                          getfem::fem_descriptor("FEM_QK(2, 0)"));
   //}
   size_type nb_dof_rhs = mf_rhs.nb_dof();
   plain_vector F(nb_dof_rhs * 3);
   plain_vector M(nb_dof_rhs * 2);
        //
   cout << "Number of dofs for ut: " << mf_ut.nb_dof() << endl;
        cout << "Number of dofs for u3: " << mf_u3.nb_dof() << endl;
   cout << "Number of dofs for theta: " << mf_theta.nb_dof() << endl;

   // define a brick
   getfem::mdbrick_abstract<> *ELAS, *SIMPLE(0);
   // plate
   getfem::mdbrick_isotropic_linearized_plate<>
   ELAS1(mim, mim_subint, mf_ut, mf_u3, mf_theta, lambda, mu, h);
   if(mitc) ELAS1.set_mitc();
   ELAS = &ELAS1;
   getfem::mdbrick_plate_source_term<>VOL_F(*ELAS,mf_rhs,F,M);
   getfem::mdbrick_plate_clamped_support<> SIMPLE1
   (VOL_F, 1, 0, getfem::AUGMENTED_CONSTRAINTS);
   SIMPLE = &SIMPLE1;
   getfem::mdbrick_plate_closing<> final_model(*SIMPLE, 0, 1);    
        cout << "Total number of variables : " << final_model.nb_dof() << endl;
        

   return EXIT_SUCCESS;
}