|
From: | Andriy Andreykiv |
Subject: | [Getfem-users] new_interpolated_fem problem |
Date: | Wed, 28 Feb 2007 14:20:05 +0100 |
User-agent: | KMail/1.9.5 |
Good day gentlemen, I'm having problems with interpolating a finite element method from one mesh to the other (getfem::new_interpolated_fem). What I want to do is the following, below there is a code that first defines two quad elements (in 2D space) and a line elemnet that crosses it arbitrarily in the middle (there is a picture bellow). I want to calculate a mass matrix (used for fictitious domain method) which is evaluated on that segment element, but in fact is calculated as a product of shape functions from the quads and the segment element (and the integration is over the segment): M=\int_{segment} Base(the setment) X Base(the quads) dx I'm getting an error: interpolated_fem::update from context : convex_index = [0] ============================================ | An error has been detected !!! | ============================================ Error in ./gmm_opt.h, line 76 : non invertible matrix I would really appreciate the help. Best regards, Andriy -- //++++++++++++++++++++++THE CODE++++++++++++++++++++++++ #include <getfem_model_solvers.h> #include <gmm.h> #include <getfem_interpolated_fem.h> /* some Getfem++ types that we will be using */ using bgeot::base_small_vector; /* special class for small (dim<16) vectors */ using bgeot::base_node; /* geometrical nodes(derived from base_small_vector)*/ using bgeot::scalar_type; /* = double */ using bgeot::size_type; /* = unsigned long */ using bgeot::base_matrix; /* small dense matrix. */ /* definition of some matrix/vector types. * default types of getfem_model_solvers.h */ typedef getfem::modeling_standard_sparse_vector sparse_vector; typedef getfem::modeling_standard_sparse_matrix sparse_matrix; typedef getfem::modeling_standard_plain_vector plain_vector; int main(int argc, char *argv[]) { try { // Building a mesh of two quads and a segment that crosses them in the middle // ___________________________________ // | | | // | | | // | | | // | +-----------------------------+ | // | | | // | | | // | | | // |_________________|__________________| // getfem::mesh quad_mesh; getfem::mesh line_mesh; std::vector<bgeot::size_type> ind1(4); std::vector<bgeot::size_type> ind2(4); std::vector<bgeot::size_type> ind3(2); ind1[0]=quad_mesh.add_point(bgeot::base_node(0.0,0.0)); ind1[1]=quad_mesh.add_point(bgeot::base_node(0.0,1.0)); ind1[2]=quad_mesh.add_point(bgeot::base_node(1.0,1.0)); ind1[3]=quad_mesh.add_point(bgeot::base_node(1.0,0.0)); ind2[0]=ind1[2]; ind2[1]=ind1[3]; ind2[2]=quad_mesh.add_point(bgeot::base_node(2.0,0.0)); ind2[3]=quad_mesh.add_point(bgeot::base_node(2.0,1.0)); ind3[0]=line_mesh.add_point(bgeot::base_node(0.5,0.5)); ind3[1]=line_mesh.add_point(bgeot::base_node(1.5,0.5)); quad_mesh.add_parallelepiped(2,ind1.begin()); quad_mesh.add_parallelepiped(2,ind2.begin()); line_mesh.add_convex(bgeot::simplex_geotrans(1,1),ind3.begin());
// mesh_fem for the two quad elements getfem::mesh_fem mq(quad_mesh); mq.set_finite_element(mq.linked_mesh().convex_index(),getfem::fem_descriptor("FEM_QK(2,1)")); mq.set_qdim(1); // mesh_fem and integration method for the segment element getfem::mesh_fem ml(line_mesh); ml.set_finite_element(getfem::fem_descriptor("FEM_PK(1,1)")); ml.set_qdim(1); getfem::mesh_im mim(line_mesh); getfem::pintegration_method ppi=getfem::int_method_descriptor("IM_GAUSS1D(3)"); mim.set_integration_method(ml.linked_mesh().convex_index(),ppi); // CREATION OF THE INTERPOLATION MEHS_FEM, THAT IS SUPPOSED TO INTERPOLATE // THE QUAD'S SHAPE FUNCTIONS ON THE SEGMENT getfem::mesh_fem mq_interpolate(mq.linked_mesh()); getfem::pfem ifem=getfem::new_interpolated_fem(mq,mim); //THIS DOESN'T // WORK, DON'T KNOW WHY !!!!!! mq_interpolate.set_finite_element(ifem); // mass matrix on the segment created from the shape function of the quads and the segment sparse_matrix M(6,2); getfem::asm_mass_matrix(M,mim,mq_interpolate,ml); // std::cout<<" M="<<M<<std::endl; // } DAL_STANDARD_CATCH_ERROR; return 0; } |
[Prev in Thread] | Current Thread | [Next in Thread] |