On Tue, Jul 13, 2010 at 02:59:35PM +0200, Christian Fischer wrote: > Dear All, > > I recently thought about using the following structural elements: > > * 1D: Euler-Bernoulli or Timoshenko beam > * 2D: axially symmetric continuum elements > * 1D: axially symmetric flat plate (like an axially symmetric beam) > > Before I start working on any of these elements. Maybe there is somebody who > already implemented any of these. (Like the truss element in a recent email.) > Or somebody has a few hints for me e.g. how to implement axial symmetry in > general. > > Best wishes > Christian Fischer > Dear Christian
I am glad for your posting. I will gladly discuss implementation of the above elements as any such discussion greatly enhances documentation of GetFEM internals. Recently I have implemented 2D frame element. I have attached some source code (the implementation of frame brick is not complete in the sense that I do not account for element stiffness). It is sparse of comments but I can serve any info. >From my understanding of GetFEM (Yves can correct me for sure :) >implementation of any of your above problems does not really require adding new elements to GetFEM. Everything can be implemented on the assembly level. For beam elements there already Hermite element (in case we agree to skip axsial displacement - if not it should be possible to handle it via separate mesh_fem). Of course having true structural elements would be nice to. I have not implemented axisymmetric problem but If I recall correctly in some post I have seen comment that functions like Lame coefficients can be implemented via mesh_fem_global_function. So again it should suffice to write right assembly routines. That is my five cents. Regards, Roman -- Roman Putanowicz, PhD < [email protected] > Institute for Computational Civil Engng (L-5) Dept. of Civil Engng, Cracow Univ. of Technology www.l5.pk.edu.pl, tel. +48 12 628 2569, fax 2034
/* * Author: Roman Putanowicz <[email protected]> * Created: Thu May 06 10:22:55 2010 * Modified by: */ /** * @file Implementation of frame element * */ #include <stdlib.h> // for EXIT_SUCCESS #include <stdio.h> #include <iostream> #include <vector> #include <map> #include <string> #include "getfem/getfem_model_solvers.h" #include "getfem/getfem_models.h" #include "getfem/getfem_export.h" #include "getfem/getfem_regular_meshes.h" #include "getfem/getfem_mesh_slice.h" #include "getfem/bgeot_comma_init.h" #include "getfem/getfem_mesh_fem_sum.h" #include "getfem/getfem_mesh_fem_global_function.h" #include "gmm/gmm.h" #include "getfem/getfem_fem.h" #include <string> namespace femdk { /* ******************************************************************** */ /* Frame element on the segment */ /* ******************************************************************** */ struct frame_segment__ : public getfem::fem<getfem::base_poly> { virtual void mat_trans(getfem::base_matrix &M, const getfem::base_matrix &G, bgeot::pgeometric_trans pgt) const; frame_segment__(void); }; void frame_segment__::mat_trans(getfem::base_matrix &M, const getfem::base_matrix &G, bgeot::pgeometric_trans pgt) const { static bgeot::pgeotrans_precomp pgp; static bgeot::pgeometric_trans pgt_stored = 0; static getfem::base_matrix K(1, 1); static getfem::base_vector r(1); getfem::dim_type N = getfem::dim_type(G.nrows()); gmm::resize(r, N); for (getfem::size_type i = 0; i < N; ++i) r[i] = ::exp(double(i)); if (pgt != pgt_stored) { pgt_stored = pgt; pgp = bgeot::geotrans_precomp(pgt, node_tab(0), 0); gmm::resize(K, N, 1); } gmm::copy(gmm::identity_matrix(), M); // gradient at point 0 gmm::mult(G, pgp->grad(1), K); if (N == 1) M(3, 3) = K(0,0); else M(3, 3) = gmm::mat_euclidean_norm(K) * gmm::sgn(gmm::vect_sp(gmm::mat_col(K, 0), r)); M(2,2) = K[0]/M(3,3); M(0,0) = K[0]/M(3,3); M(0,2) = K[1]/M(3,3); M(2,0) = K[1]/M(3,3); // gradient at point 1 if (!(pgt->is_linear())) { gmm::mult(G, pgp->grad(3), K); } if (N == 1) M(5, 5) = K(0,0); else M(5, 5) = gmm::mat_euclidean_norm(K) * gmm::sgn(gmm::vect_sp(gmm::mat_col(K, 0), r)); M(4,4) = K[0]/M(5,5); M(1,1) = K[0]/M(5,5); M(1,4) = K[1]/M(5,5); M(4,1) = K[1]/M(5,5); } frame_segment__::frame_segment__(void) { getfem::base_node pt(1); cvr = bgeot::simplex_of_reference(1); dim_ = cvr->structure()->dim(); init_cvs_node(); es_degree = 3; is_pol = true; is_lag = false; is_equiv = false; ntarget_dim = 2; base_.resize(12); for (getfem::size_type j = 0; j < getfem::size_type(12); ++j) base_[j] = bgeot::null_poly(1); pt[0] = 0.0; add_node(getfem::to_coord_dof(getfem::lagrange_dof(1), 0), pt); getfem::read_poly(base_[0], 1, "(1 - x)"); pt[0] = 1.0; add_node(getfem::to_coord_dof(getfem::lagrange_dof(1), 0), pt); getfem::read_poly(base_[1], 1, "x"); pt[0] = 0.0; add_node(getfem::to_coord_dof(getfem::lagrange_dof(1), 1), pt); getfem::read_poly(base_[8], 1, "(1 - x)^2*(2*x + 1)"); pt[0] = 0.0; add_node(getfem::to_coord_dof(getfem::derivative_dof(1, 0), 1), pt, dal::bit_vector()); getfem::read_poly(base_[9], 1, "x*(x - 1)*(x - 1)"); pt[0] = 1.0; add_node(getfem::to_coord_dof(getfem::lagrange_dof(1), 1), pt); getfem::read_poly(base_[10], 1, "x*x*(3 - 2*x)"); pt[0] = 1.0; add_node(getfem::to_coord_dof(getfem::derivative_dof(1, 0), 1), pt, dal::bit_vector()); getfem::read_poly(base_[11], 1, "x*x*(x - 1)"); } static getfem::pfem Frame_fem(getfem::fem_param_list ¶ms, std::vector<dal::pstatic_stored_object> &dependencies) { GMM_ASSERT1(params.size() == 1, "Bad number of parameters : " << params.size() << " should be 1."); GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters"); int d = int(::floor(params[0].num() + 0.01)); getfem::virtual_fem *p = 0; switch(d) { case 1 : p = new frame_segment__; break; case 2 : p = new frame_segment__; break; default : GMM_ASSERT1(false, "Sorry, Frame element in dimension " << d << " not available"); } dependencies.push_back(p->ref_convex(0)); dependencies.push_back(p->node_tab(0)); return p; } void use_frame_elem() { dal::singleton<getfem::fem_naming_system>::instance().add_suffix("FRAME", Frame_fem); } } /* namespac femdk */
/* * Author: Roman Putanowicz <[email protected]> * Created: Thu May 06 10:22:55 2010 * Modified by: */ /** * @file Bending of a simple beam * */ # define GMM_TRACES_LEVEL 0 #include "structural_fem_config.h" extern "C" { #include "simple_beam_options.h" } #include <stdlib.h> // for EXIT_SUCCESS #include <stdio.h> #include <iostream> #include <vector> #include <map> #include <string> #include "getfem/dal_bit_vector.h" #include "getfem/getfem_model_solvers.h" #include "getfem/getfem_models.h" #include "getfem/getfem_modeling.h" #include "getfem/getfem_export.h" #include "getfem/getfem_regular_meshes.h" #include "getfem/getfem_mesh_slice.h" #include "getfem/bgeot_comma_init.h" #include "getfem/getfem_mesh_fem_sum.h" #include "getfem/getfem_mesh_fem_global_function.h" #include "gmm/gmm.h" #include "gmm/gmm_inoutput.h" #include "getfem/getfem_fem.h" #include "getfem/getfem_assembling_tensors.h" #include <string> #include "femdk/femdk_meshing.hxx" namespace femdk { extern void use_frame_elem(); } namespace getfem { struct frame_brick : public virtual_brick { void asm_real_tangent_terms(const model &md, getfem::size_type ib, const model::varnamelist &vl, const model::varnamelist &dl, const model::mimlist &mims, model::real_matlist &matl, model::real_veclist &vecl, model::real_veclist &vecl_sym, size_type region, model::build_version version) const { GMM_ASSERT1(mims.size() == 1, "Frame brick need a single mesh_im"); GMM_ASSERT1(vl.size() == 1, "Frame brick need a single variable"); GMM_ASSERT1(dl.size() == 1, "Wrong number of data for frame brick, " << dl.size() << " should be 1 (vector)."); GMM_ASSERT1(matl.size() == 1, "Wrong number of terms for frame brick "); const model_real_plain_vector &u = md.real_variable(vl[0]); const mesh_fem &mf_u = *(md.pmesh_fem_of_variable(vl[0])); const mesh_fem *mf_params = md.pmesh_fem_of_variable(dl[0]); const model_real_plain_vector ¶ms = md.real_variable(dl[0]); const mesh_im &mim = *mims[0]; mesh_region rg(region); mf_u.linked_mesh().intersect_with_mpi_region(rg); if (version & model::BUILD_MATRIX) { gmm::clear(matl[0]); GMM_TRACE2("Frame stiffness matrix assembly"); getfem::generic_assembly assem; assem.push_mi(mim); assem.push_mf(mf_u); assem.push_mat(matl[0]); assem.set("t=comp(vHess(#1).vHess(#1));" "g=comp(vGrad(#1).vGrad(#1));" "Z=sym(t(:,2,k,k,:,2,j,j) + g(:,1,i,:,1,i));" "M(#1,#1)+=Z;" ); assem.assembly(); } } /* end of assem_real_tangent_terms */ frame_brick() { set_flags("Frame brick", true /* is linear*/, true /* is symmetric */, true /* is coercive */, true /* is real */, false /* is complex */); } }; size_type add_frame_brick (model &md, const mesh_im &mim, const std::string &varname, const std::string &dataname, size_type region) { pbrick pbr = new frame_brick(); model::termlist tl; tl.push_back(model::term_description(varname, varname, true)); model::varnamelist dl(1, dataname); model::varnamelist vl(1, varname); return md.add_brick(pbr, vl, dl, tl, model::mimlist(1,&mim), region); } } /* end of namespace getfem */ int main_body(int argc, char** argv); /** main function * * This function should contain command line parsing and initialization * code. The real work should be done by calling main_body(). */ int main(int argc, char **argv) { int result = EXIT_SUCCESS; result = main_body(argc, argv); return result; } /** The function where the main work is done after configuration * and command line parsing. */ int main_body(int argc, char** argv) { try { { int optct = optionProcess( &simple_beamOptions, argc, argv ); argc -= optct; argv += optct; } /* Benchmark data */ double L = atof(OPT_ARG(L)); double N = atoi(OPT_ARG(N)); femdk::use_frame_elem(); getfem::mesh mesh; getfem::model model; femdk::meshing_params params; params.N = N ; femdk::make_segment_mesh(mesh, bgeot::base_node(0.0,0.0,0.0), bgeot::base_node(L, 0.0,0.0), params); /* Build regions corresponding to the end points */ getfem::mesh_region leftEnd = mesh.region(1); getfem::mesh_region rightEnd = mesh.region(2); /* Find number of points in the mesh */ getfem::size_type nnodes = mesh.nb_points(); getfem::size_type leftNode = 0; getfem::size_type rightNode = nnodes-1; getfem::size_type leftConvex = mesh.first_convex_of_point(leftNode); getfem::size_type rightConvex = mesh.first_convex_of_point(rightNode); getfem::size_type leftFace = 1; getfem::size_type rightFace = 0; leftEnd.add(leftConvex, leftFace); rightEnd.add(rightConvex, rightFace); getfem::mesh_fem femFrame(mesh); femFrame.set_qdim(2); femFrame.set_finite_element(getfem::fem_descriptor("FEM_FRAME(2)")); getfem::mesh_im mim(mesh); getfem::pintegration_method pim = getfem::int_method_descriptor("IM_GAUSS1D(2)"); mim.set_integration_method(mesh.convex_index(), pim); getfem::size_type Ndof = femFrame.nb_dof(); std::cout << "NUMBER OF DEGREES OF FREEDOM " << Ndof << "\n"; typedef gmm::wsvector<getfem::scalar_type> model_real_sparse_vector; typedef std::vector<getfem::scalar_type> model_real_plain_vector; typedef gmm::col_matrix<model_real_sparse_vector> model_real_sparse_matrix; typedef gmm::rsvector<getfem::scalar_type> sparse_vector_type; typedef gmm::row_matrix<sparse_vector_type> sparse_matrix_type; typedef gmm::col_matrix<sparse_vector_type> col_sparse_matrix_type; typedef std::vector<getfem::scalar_type> plain_vector; model.add_fem_variable("u", femFrame); plain_vector frame_load(2, 0.0); frame_load[1] = 12.0; model.add_initialized_fixed_size_data("load", frame_load); getfem::add_source_term_brick(model, mim, "u", "load"); plain_vector frame_prop(3, 0.0); model.add_initialized_fixed_size_data("EIA", frame_prop); getfem::add_frame_brick(model,mim, "u", "EIA", 0); plain_vector bc(2, 0.0); plain_vector bc1(2, 5.0); model.add_initialized_fixed_size_data("bc", bc); model.add_initialized_fixed_size_data("bc1", bc1); // getfem::add_source_term_brick(model, mim, "u", "bc1",1); // getfem::add_source_term_brick(model, mim, "u", "bc1",2); getfem::add_Dirichlet_condition_with_penalization(model, mim, "u", 1e9, 1, "bc"); getfem::add_Dirichlet_condition_with_penalization(model, mim, "u", 1e9, 2, "bc"); plain_vector U; gmm::resize(U, Ndof); gmm::clear(U); gmm::iteration iter(1.0e-5, 0, 55); iter.init(); getfem::standard_solve(model, iter, true); gmm::copy(model.real_variable("u"), U); double dx = 1.0/params.N; int j; for (unsigned int i=1,j=0; i<U.size(); i+=3, j++) { std::cout << j*dx << " " << U[i] << "\n"; } return EXIT_SUCCESS; } GMM_STANDARD_CATCH_ERROR; return EXIT_FAILURE; }
_______________________________________________ Getfem-users mailing list [email protected] https://mail.gna.org/listinfo/getfem-users
