Dear Jano,
IM_GAUSS1D(3) represent a (small) sub-integration. Did you try with IM_GAUSS1D(5) ? Yves. Jano <[email protected]> a écrit : > Hi, > > I just started using the getfem++ lib, made some tests but I'm getting > some weird results. > I wanted to test the equation: > d4u/dx4 = f > u = variable > f = constant = 2 > b.c. > > u(0) = u(1) = 0 > du(0) = du(1) = 0 > > using FEM_HERMITE(1), simplex_geotrans(1,1) and IM_GAUSS1D(3), > > for what I know it is a well conditioned problem and should give the > exact solution. > It compiles, run, but give a huge condition number and nothing even > near the exact solution. > > thanks if anyone can help, > Jano > > ps: sorry if the email is duplicated > > > Here is my code if someone can find what i did wrong: > > #include <stdlib.h> > #include <string> > #include <getfem/getfem_mesh.h> > #include <getfem/getfem_mesh_fem.h> > #include <getfem/getfem_import.h> > #include <getfem/getfem_regular_meshes.h> > #include <getfem/getfem_mesh_im.h> > #include <getfem/getfem_models.h> > #include <gmm/gmm.h> > #include <getfem/getfem_export.h> > #include <getfem/getfem_model_solvers.h> > #include <getfem/getfem_fourth_order.h> > > bgeot::scalar_type condicaoDirichlet(const bgeot::base_node &x) { > return 0.0; > } > > bgeot::scalar_type DcondicaoDirichlet(const bgeot::base_node &x) { > return 0.0; > } > > bgeot::scalar_type D = 1.0; > bgeot::scalar_type f = -2.0; > > > /* > * > */ > int main(int argc, char** argv) { > > //FEM_HERMITE(n)-Elemento de Hermite com n dimensões. > std::string femName = "FEM_HERMITE(1)"; > //Flag para a regiao de contorno > int DIRICHLET_BOUNDARY_NUM = 0; > //Malha > getfem::mesh malha; > //getfem::import_mesh(string FILENAME, string FORMAT, mesh m) > //Geracao e refinamento da malha > std::vector<getfem::size_type> nsubdiv(1); > nsubdiv[0] = 50; > //nsubdiv[1] = 100; > getfem::regular_unit_mesh(malha, nsubdiv, bgeot::simplex_geotrans(1,1)); > //mesh_fem representando U e f > getfem::mesh_fem meshFem(malha); > getfem::mesh_fem meshFemRightHandSide(malha); > //Define o tipo de elemento finito > meshFem.set_finite_element(getfem::fem_descriptor(femName)); > meshFemRightHandSide.set_finite_element(getfem::fem_descriptor(femName)); > //IM_NC(n,k)-Integracao numerica com Newton-Cotes > getfem::mesh_im meshIm(malha); > > meshIm.set_integration_method(getfem::int_method_descriptor("IM_GAUSS1D(3)")); > //Instanciacao do modelo > getfem::model modelo; > //Atribui U e o metodo de integracao > modelo.add_fem_variable("u", meshFem); > //getfem::add_Laplacian_brick(modelo, meshIm, "u"); > modelo.add_initialized_scalar_data("D", D); > getfem::add_bilaplacian_brick(modelo, meshIm, "u", "D"); > modelo.add_initialized_scalar_data("f", f); > getfem::add_source_term_brick(modelo, meshIm, "u", "f"); > //Define a regiao de contorno com a flag DIRICHLET_BOUNDARY_NUM > getfem::mesh_region border_faces; > getfem::outer_faces_of_mesh(malha, border_faces); > for (getfem::mr_visitor i(border_faces); !i.finished(); ++i) { > malha.region(DIRICHLET_BOUNDARY_NUM).add(i.cv(), i.f()); > } > //Condicoes de contorno de Dirichlet > std::vector<bgeot::scalar_type> F(meshFemRightHandSide.nb_dof()); > std::vector<bgeot::scalar_type> F2(meshFemRightHandSide.nb_dof()); > getfem::interpolation_function(meshFemRightHandSide, F, > condicaoDirichlet); > getfem::interpolation_function(meshFemRightHandSide, F2, > DcondicaoDirichlet); > modelo.add_initialized_fem_data("DirichletData", > meshFemRightHandSide, F); > modelo.add_initialized_fem_data("derivadaDirichletData", > meshFemRightHandSide, F2); > getfem::add_Dirichlet_condition_with_multipliers( > modelo, meshIm, "u", meshFem, DIRICHLET_BOUNDARY_NUM, > "DirichletData"); > getfem::add_normal_derivative_Dirichlet_condition_with_multipliers( > modelo, meshIm, "u", meshFem, DIRICHLET_BOUNDARY_NUM, > "derivadaDirichletData"); > > > //Define o residuo maximo e um iterator > bgeot::scalar_type maxResidualIterativeSolver = pow(10, -5); > gmm::iteration iter(maxResidualIterativeSolver, 1, 40000); > //Solver > getfem::standard_solve(modelo, iter); > //Guarda a solucao em u > std::vector<bgeot::scalar_type> u = modelo.real_variable("u"); > > std::fstream solucao("solucao.txt",std::ios::out); > for (unsigned i=0; i < gmm::vect_size(u); ++(++i)) > solucao << u[i] << "\t" << u[i+1] << "\n"; > > // when the 2nd arg is true, the mesh is saved with the |mf| > meshFem.write_to_file("solucao.mf", true); > > > //Exporta > getfem::dx_export exp_dx("solucao.dx"); > exp_dx.exporting(malha); > exp_dx.write_point_data(meshFem, u, "deslocamento"); > > getfem::vtk_export exp_vtk("solucao.vtk"); > exp_vtk.exporting(malha); > exp_vtk.write_point_data(meshFem, u, "deslocamento"); > > return (EXIT_SUCCESS); > } > > _______________________________________________ > 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
