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
_______________________________________________
Getfem-users mailing list
[email protected]
https://mail.gna.org/listinfo/getfem-users
