#include "getfem/getfem_mesh.h"
#include "getfem/getfem_mesh_fem.h"
#include "getfem/getfem_import.h"
#include "getfem/getfem_generic_assembly.h"
#include "gmm/gmm.h"
#include "getfem/getfem_model_solvers.h"
#include "getfem/getfem_mesher.h"
#include "getfem/getfem_fem.h"
using std::endl; using std::cout;
using bgeot::base_small_vector;
typedef getfem::model_real_plain_vector plain_vector;
int main(){
getfem::mesh m;
getfem::import_mesh("gmsh:./Simple.msh",m);
getfem::mesh_region border_faces;
getfem::outer_faces_of_mesh(m, border_faces);
getfem::mesh_region fb1 = getfem::select_faces_of_normal(m, border_faces,
base_small_vector( 1., 0.), 0.01);
getfem::mesh_region fb2 = getfem::select_faces_of_normal(m, border_faces,
base_small_vector( -1., 0.), 0.01);
m.region(1) = fb1; m.region(2) = fb2;
// Set finite element methods
getfem::pfem pf = getfem::fem_descriptor("FEM_QK(2,1)");
getfem::mesh_fem mfu(m); // Finite element for the elastic displacement
mfu.set_qdim(2);
mfu.set_finite_element(pf);
getfem::mesh_im mim(m); // Integration method
mim.set_integration_method(bgeot::dim_type(2));
getfem::model md;
md.add_fem_variable("u", mfu); // Displacement of the structure
// Membrane elastic deformation
double E = 21E6;
double nu = 0.3;
double F = 100E2;
double clambda = E*nu/((1+nu)*(1-2*nu)); // First Lame coefficient (N/cm^2)
double cmu = E/(2*(1+nu)); // Second Lame coefficient (N/cm^2)
double clambdastar = 2*clambda*cmu/(clambda+2*cmu); // Lame coefficient
md.add_initialized_scalar_data("cmu", cmu);
md.add_initialized_scalar_data("clambdastar", clambdastar);
getfem::add_isotropic_linearized_elasticity_brick
(md, mim, "u", "clambdastar", "cmu");
getfem::add_Dirichlet_condition_with_multipliers
(md, mim, "u", bgeot::dim_type(2-1), 2);
md.add_initialized_fixed_size_data("Fdata", base_small_vector(F,0.));
getfem::add_source_term_brick(md, mim, "u", "Fdata", 1);
gmm::iteration iter(1E-9, 1, 100);
getfem::standard_solve(md, iter);
//
// Solution export
plain_vector U(mfu.nb_dof());
gmm::copy(md.real_variable("u"), U);
bgeot::base_vector coeff;
getfem::slice_vector_on_basic_dof_of_element(mfu,U,1,coeff);
std::vector<std::string> v = { "Point a ", "Point b ", "Point c ", "Point d "};
cout << "Displacements for element 1: \n";
for (int i=0; i < coeff.size()/2; i++)
cout << v[i] << "u1: " << coeff[2*i] << "; u2: " << coeff[2*i+1] << endl;
return 0;
}