[Getfem-users] new_interpolated_fem problem

[Andriy Andreykiv]

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;

}