Hi Tom,
if you have compiled GetFEM++ from source maybe you can try replacing:
asm_real_or_complex_1_param
(*B, mim, mf_mult, *mf_H, *H, rg, (mf_u.get_qdim() == 1) ?
"F=data(#2);"
"M(#1,#3)+=comp(Base(#1).Base(#3).Base(#2)(:,:,i).F(i))"
: "F=data(qdim(#1),qdim(#1),#2);"
"M(#1,#3)+=comp(vBase(#1).vBase(#3).Base(#2)(:,i,:,j,k).F(i,j,k));", &mf_u);
in Dirichlet_condition_brick::asm_real_tangent_terms in file
getfem_models.cc, with:
asm_real_or_complex_1_param
(*B, mim, mf_mult, *mf_H, *H, rg, (mf_u.get_qdim() == 1) ?
"F=data(#2);"
"M(#1,#3)+=comp(Base(#1).Base(#3).Base(#2))(:,:,i).F(i)"
: "F=data(qdim(#1),qdim(#1),#2);"
"M(#1,#3)+=comp(vBase(#1).vBase(#3).Base(#2))(:,i,:,j,k).F(i,j,k);", &mf_u);
and recompile.
BR
Kostas
On Thu, Jan 8, 2015 at 3:14 PM, Tom Haeck <[email protected]>
wrote:
> 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
>
>
_______________________________________________
Getfem-users mailing list
[email protected]
https://mail.gna.org/listinfo/getfem-users
