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
address@hidden
https://mail.gna.org/listinfo/getfem-users