Hi!
I'm experimenting with integration on faces using the high-level
assembly. I was not expecting that the following two expressions would
give different results when integrating on a single face in a
FEM_QK(2,1) mesh:
Case 1: "Test_u*Grad_u.Normal"
Case 2: "Test_u*Normal.Grad_u"
For all the details, please see the attached test program which shows
that these expressions give different results in the case of integrating
along face 0 in convex 0, in a 2x2 FEM_QK(2,1) mesh. I expect Case 1 to
be the correct behaviour, because its results coincide with the result
of selecting the x-component of Grad_u explicitly (since the outward
unit normal for the face in question is [1,0]):
Case 3: "Test_u*Grad_u(1)"
So unless I have misunderstood, something wrong happens in Case 2. I'm
using a very recent SVN-version of Getfem++ and the latest GCC 5 snapshot.
Best regards,
Torquil Sørensen
#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);
}
_______________________________________________
Getfem-users mailing list
[email protected]
https://mail.gna.org/listinfo/getfem-users
