#include "getfem/getfem_assembling.h"
#include "getfem/getfem_generic_assembly.h"
#include "getfem/getfem_mesh.h"
#include "getfem/getfem_mesh_fem.h"
#include "getfem/getfem_mesh_im.h"
#include "getfem/getfem_mesh_region.h"
#include "getfem/getfem_regular_meshes.h"

int main()
{
	getfem::mesh mesh;
	std::vector<bgeot::size_type> ref(2); ref[0] = 2; ref[1] = 2;

	getfem::regular_unit_mesh(mesh, ref, bgeot::parallelepiped_geotrans(2, 1));

	getfem::mesh_fem mf(mesh);

	mf.set_finite_element(getfem::fem_descriptor("FEM_QK(2,1)"));

	getfem::mesh_im mi(mesh);
	mi.set_integration_method(getfem::int_method_descriptor("IM_GAUSS_PARALLELEPIPED(2,4)"));

	// High-level assembly of integral of n.(v*grad_u) on face (0,0):

	getfem::mesh_region reg_face_00;
	reg_face_00.add(0,0);

	{	// Case 1
		getfem::model_real_sparse_matrix mat(mf.nb_dof(), mf.nb_dof());
		getfem::ga_workspace workspace;
		getfem::base_vector u(mf.nb_dof());
		workspace.add_fem_variable("u", mf, gmm::sub_interval(0, mf.nb_dof()), u);
		workspace.add_expression("Test_u*Grad_u.Normal", mi, reg_face_00);
		workspace.set_assembled_matrix(mat);
		workspace.assembly(2);
		std::cout << "mat (case 1) = " << mat << std::endl;
	}

	{	// Case 2
		getfem::model_real_sparse_matrix mat(mf.nb_dof(), mf.nb_dof());
		getfem::ga_workspace workspace;
		getfem::base_vector u(mf.nb_dof());
		workspace.add_fem_variable("u", mf, gmm::sub_interval(0, mf.nb_dof()), u);
		workspace.add_expression("Test_u*Normal.Grad_u", mi, reg_face_00);
		workspace.set_assembled_matrix(mat);
		workspace.assembly(2);
		std::cout << "mat (case 2) = " << mat << std::endl;
	}

	{	// Case 3 (selecting the x-component of Grad_u explicitly, since the unit normal vector should be [1,0])
		getfem::model_real_sparse_matrix mat(mf.nb_dof(), mf.nb_dof());
		getfem::ga_workspace workspace;
		getfem::base_vector u(mf.nb_dof());
		workspace.add_fem_variable("u", mf, gmm::sub_interval(0, mf.nb_dof()), u);
		workspace.add_expression("Test_u*Grad_u(1)", mi, reg_face_00);
		workspace.set_assembled_matrix(mat);
		workspace.assembly(2);
		std::cout << "mat (case 3) = " << mat << std::endl;
	}

	return(0);
}