Re: [Getfem-users] plate elements

[Yves Renard]

Dear Umut,

In order to compute the stiffness matrix, you have to use a model_state 
variable MS

 MS.adapt_sizes(final_model);
 final_model.compute_tangent_matrix(MS);
 final_model.compute_residual(MS);

and eventually

MS.compute_reduced_system();

The stiffness matrix is then given by 

MS.tangent_matrix()

or

MS.reduced_tangent_matrix()



If you do not really use the brick system, you can also directly use the 
assembly procedures given in getfem/getfem_linearized_plates.h


The mass matrix has not been implemented yet (not very complex but we did not 
study dynamic plate models). 


Yves.



On samedi 5 septembre 2009, Umut Tabak wrote:
> 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;
> }
>
>
>
> _______________________________________________
> Getfem-users mailing list
> address@hidden
> https://mail.gna.org/listinfo/getfem-users



-- 

  Yves Renard (address@hidden)       tel : (33) 04.72.43.87.08
  Pole de Mathematiques, INSA-Lyon             fax : (33) 04.72.43.85.29
  20, rue Albert Einstein
  69621 Villeurbanne Cedex, FRANCE
  http://math.univ-lyon1.fr/~renard

---------