Dear all,
I was trying to add the pointwise constraints for a pure Neumann problem. In elastostatic.cc example, I eliminated the Dirichelet bc's and added pointwise constraints, see below
/***************************************************************/
/* top and bottom are selected as Neumann bc's*/
void elastostatic_problem::select_boundaries(void) {
size_type N = mesh.dim();
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 (gmm::abs(un[0] - 0.0) < 0.1) { // new Neumann face
mesh.region(NEUMANN_BOUNDARY_NUM).add(
i.cv(), i.f());
}
}
}
/***************************************************************/
/***************************************************************/
/* add pointwise constraint
std::vector<scalar_type> cpoints(4);
cpoints[0] = 0.0; cpoints[1] = 0.0;
cpoints[2] = 0.0; cpoints[3] = 0.0;
std::vector<scalar_type> cunitv(4);
cunitv[0] = 1.0; cunitv[1] = 0.0;
cunitv[2] = 0.0; cunitv[3] = 1.0;
model.add_initialized_fixed_size_data("cpoints", cpoints);
model.add_initialized_fixed_size_data("cunitv",cunitv);
getfem::add_pointwise_constraints_with_multipliers(model,"u","cpoints","cunitv");
/***************************************************************/
The output is
/***************************************************************/
Number of dof for u: 51842
Total number of variables : 51844
Trace 2 in ../getfem_models.cc, line 3265: Stiffness matrix assembly for isotropic linearized elasticity
Trace 2 in ../getfem_models.cc, line 1522: Source term assembly
Trace 2 in ../getfem_models.cc, line 1657: source term assembly
Level 1 Warning in ../getfem_superlu.cc, line 216: SuperLU solve failed: info =51845condition number: 8.47747e+17
terminate called after throwing an instance of 'gmm::gmm_error'
what(): Error in ../elastostatic.cc, line 542 int main(int, char**):
Solve has failed
/***************************************************************/
It seems that the rigid displacements are not correctly enforced. Could anyone tell me where I was wrong? Thanks.
Regards,
Wen