[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Getfem-users] Basic example problem

From: Renard Yves
Subject: Re: [Getfem-users] Basic example problem
Date: Tue, 24 Mar 2009 20:55:30 +0100
User-agent: Dynamic Internet Messaging Program (DIMP) H3 (1.1)

Dear Daniel,

I think this is because of your last loop. You move the point of your mesh with U, but U is the value of dofs which have their own numbering (moreover you use a P2 finite element). So you have first to establish the correspondance between the node indices and the dof indices. You can use for that the function mf.ind_dof_of_element(...) or identify the point using mf.point_of_dof(...) or use an interpolation on a P1 finite element if you want to make a generic function ...


Daniel Burkhart <address@hidden> a écrit :

Hi everyone,
I have tried to convert the tripod example (from from the
matlab interface to the C++ interface.

The attached Screenshots show input mesh (green region is TOP, boundary),
and the (very noisy) output mesh)

Any suggestion about errors in my code would be great. Maybe the last lines,
where I update the positions of the nodes are incorrect, but I don't know
the correct choice!

And here my simulation Code:

Precondition: The tripod mesh is loaded to m_mesh, and the boundaries
(TOP,BOTTOM) are initialized. This is done with a GUI.

std::cout << "MESH:" << std::endl;
    std::cout << "   Nodes: " << m_mesh.nb_points() << std::endl;
    std::cout << "   Convex: " << m_mesh.nb_convex() << std::endl;

    bool LINEAR = true;
    bool INCOMPRESSIBLE = false;
    int TOP = 0;
    int BOTTOM = 1;

    getfem::mesh_fem mfu(m_mesh, 3);  //mesh-fem supporting a 3D-vector
    getfem::mesh_fem mfd(m_mesh, 1);  //scalar mesh_fem, for data fields.

    //the mesh_im stores the integration methods for each tetrahedron
    getfem::mesh_im mim(m_mesh);


    //we choose a P2 fem for the main unknown

    //the material is homogeneous, hence we use a P0 fem for the data

    //boundary stuff
    //ensure top is boundary 0, bottom is boundary 1

    //Lamé coefficient
    double E = 1000;
    double Nu = 0.3;
    //set the Lame coefficients
    double lambda = E*Nu/((1+Nu)*(1-2*Nu));
    double mu = E/(2*(1+Nu));

    //create a meshfem for the pressure field (used if incompressible  = 0)
    getfem::mesh_fem mfp(m_mesh);


    //mesh stats:
    cout << "Number of dofs for u:" << mfu.nb_dof() << endl;
    cout << "Number of dofs for p:" << mfp.nb_dof() << endl;

    //bricks stuff

    // the linearized elasticity , for small displacements
    getfem::mdbrick_abstract<> *b0=0;
    getfem::mdbrick_abstract<> *b1=0;
    getfem::mdbrick_abstract<> *b2=0;
    getfem::mdbrick_abstract<> *b3=0;

    //getfem::mdbrick_isotropic_linearized_elasticity<> *b0 = 0;
    //getfem::mdbrick_linear_incomp<> *b1 = 0;

        b0 = new getfem::mdbrick_isotropic_linearized_elasticity<>(mim, mfu,
lambda, mu);
            b1 = new getfem::mdbrick_linear_incomp<>(*b0,mfp);
        bgeot::base_vector p(2); p[0] = 0.02; p[1] = 0.5;

            getfem::Mooney_Rivlin_hyperelastic_law HypLaw;
            b0 = new getfem::mdbrick_nonlinear_elasticity<>(HypLaw, mim,
mfu, p);
            b1 = new getfem::mdbrick_nonlinear_incomp<>(*b0,mfp);
            getfem::SaintVenant_Kirchhoff_hyperelastic_law HypLaw;
            b0 = new getfem::mdbrick_nonlinear_elasticity<>(HypLaw, mim,
mfu, p);
            b1 = b0;

    std::cout << "Init boundary conditions..." << std::endl;
    //boundary conditions:
    //set a vertical force on the top of the tripod
    bgeot::base_vector q(3); q[0] = 0.0; q[1] = -10.0; q[2]=0.0;
    b2 = new getfem::mdbrick_source_term<>(*b1,mfd, q, TOP);

    //atach to the ground
    b3 = new getfem::mdbrick_Dirichlet<> (*b2, BOTTOM);


    std::cout << "Calculation start\n";
    getfem::standard_model_state MS(*b3);

    getfem::modeling_standard_plain_vector F(mfd.nb_dof()*3);


    gmm::iteration iter(1e-10, 2, 5000);
    getfem::standard_solve(MS, *b3, iter);

    // Solution extraction
    getfem::modeling_standard_plain_vector U(mfu.nb_dof());


    std::cout << "Finished: " << iter.converged() << std::endl;

    int j=0;
    for (dal::bv_visitor i(m_mesh.points_index()); !i.finished(); ++i)
        m_mesh.points()[i][0] -= U[j++];
        m_mesh.points()[i][1] -= U[j++];
        m_mesh.points()[i][2] -= U[j++];


Thanks for your help!

reply via email to

[Prev in Thread] Current Thread [Next in Thread]