|
From: | Konstantinos Poulios |
Subject: | Re: [Getfem-users] Generalized dirichlet boundary condition for Isotropic linearised elasticity brick |
Date: | Thu, 8 Jan 2015 15:58:45 +0100 |
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 Hmodel.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 Hplain_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
_______________________________________________
Getfem-users mailing list
address@hidden
https://mail.gna.org/listinfo/getfem-users
[Prev in Thread] | Current Thread | [Next in Thread] |