[Getfem-users] Generalized dirichlet boundary condition for Isotropic linearised elasticity brick

[Tom Haeck]

Hi all,

I keep getting an error when I add a generalised dirichlet condition in the getfem::add_isotropic_linearized_elasticity_brick model.  When the matrix ‘H’ in the generalised dirichlet condition 'H*u=r' is a constant matrix, there is no problem.  However, when the matrix ‘H' is a matrix field, the problem arises.  

I created a small toy example to illustrate my problem.  A 2D disc behaves as an isotropic elastic material.  Forces act on the left hand side of the disc.  I impose a Dirichlet boundary condition on the outer points on the right hand side of the disc.  

- When I require those points to have displacement zero in all directions, i.e.  H is a fixed 2-by-2 unity matrix and r=[0,0], there is no problem.  H is declared as a plain_vector with size 4 and added to the model with the function add_initialized_fixed_size_data()

- When I require those points to have displacements zero in the direction normal to the disc, H needs to be a matrix field.  H is declared as a plain_vector with size 4*mesh_fem.nb_dof() and added to the model with the function add_initialized_fem_data().  Mesh_fem is the fem method on which the matrix field is described.  In this case, no matter what values I put in the matrix field H, I get an error:

Error in getfem_assembling_tensors.cc, line 811 : 

tensor error: wrong number of indices for the 3th argument of the reduction comp(vBase(#1).vBase(#3).Base(#2)(:,i,:,j,k).F(i,j,k)) (expected 1 indexes, got 5

The total script is as follows:

#include "getfem/getfem_export.h"

#include "getfem/getfem_model_solvers.h"

using bgeot::dim_type;

using bgeot::size_type;

using bgeot::base_node;

using bgeot::scalar_type;

typedef getfem::modeling_standard_plain_vector  plain_vector;

int main(int argc, char *argv[]) {

    

    getfem::mesh mesh;

    mesh.read_from_file("/Users/thaeck0/Software/getfem-4.3/tests/meshes/disc_2D_degree3.mesh");

    

    enum { LEFT = 0, RIGHT = 1 };

    getfem::mesh_region border_faces;

    getfem::outer_faces_of_mesh(mesh, border_faces);

    for (getfem::mr_visitor i(border_faces); !i.finished(); ++i) {

        base_node un = mesh.normal_of_face_of_convex(i.cv(), i.f());

        un /= gmm::vect_norm2(un);

        if (un[0] < 0.0) mesh.region(LEFT).add(i.cv(), i.f());

        else            mesh.region(RIGHT).add(i.cv(), i.f());

    }    

    

    scalar_type lambda = 1.0;

    scalar_type mu = 1.0;

    

    plain_vector U;

    getfem::model model;

    getfem::mesh_im mim(mesh);

    getfem::mesh_fem mf_u(mesh);

    getfem::mesh_fem mf_rhs(mesh);

    mf_u.set_qdim(dim_type(2));

    

    mim.set_integration_method(getfem::int_method_descriptor("IM_TRIANGLE(6)"));

    mf_u.set_finite_element(getfem::fem_descriptor("FEM_PK(2,1)"));

    mf_rhs.set_finite_element(getfem::fem_descriptor("FEM_PK(2,1)"));

    

    model.add_fem_variable("u", mf_u);

    model.add_initialized_scalar_data("lambda", lambda);

    model.add_initialized_scalar_data("mu", mu);

    getfem::add_isotropic_linearized_elasticity_brick(model, mim, "u", "lambda", "mu");

    bgeot::base_vector q(2); q[0] = 0.1; q[1] = 0.0;

    model.add_initialized_fixed_size_data("force", q);

    getfem::add_source_term_brick(model, mim, "u", "force",LEFT);

    // 1. Generalized Dirichlet condition with fixed H

    model.add_initialized_fixed_size_data("r_fixed",plain_vector(2, 0.0));

    plain_vector HH_fixed(4, 0.0); HH_fixed[0]=1.0; HH_fixed[3]=1.0;

    model.add_initialized_fixed_size_data("H_fixed", HH_fixed);

    getfem::add_generalized_Dirichlet_condition_with_multipliers(model, mim, "u", mf_u, RIGHT, "r_fixed", "H_fixed");

    

    // 2. Generalized Dirichlet condition with varying H

    plain_vector r_field(mf_rhs.nb_dof()*2, 0.0);

    plain_vector HH_field(mf_rhs.nb_dof()*4, 0.0);

    model.add_initialized_fem_data("r_field", mf_rhs, r_field);

    model.add_initialized_fem_data("H_field", mf_rhs, HH_field);

    getfem::add_generalized_Dirichlet_condition_with_multipliers(model, mim, "u", mf_u, RIGHT, "r_field", "H_field");

    

    gmm::iteration iter(1e-20, 1, size_type(-1));

    iter.init();

    getfem::standard_solve(model, iter);

    gmm::resize(U, mf_u.nb_dof());

    gmm::copy(model.real_variable("u"), U);

    if (!iter.converged()) GMM_ASSERT1(false, "Solve has failed. No convergence.");

    

    return 0;

    

}

Thanks!

Tom