[Getfem-users] Level Set computational time

[Gerry Rizzo]

Hi,

I'm using GetFem++ to simulate a flux in a porous media in a small 3 dimensional grid with 3x3x5 cells; I have used a mixed finite element method, with RT0 for velocity and P0 for the pressure field. The grid is obtained splitting each cell in tetrahedra.

After the problem is solved I would like to compute the mean pressure over an internal block that doesn't match the original grid; to do this I've used the Level Set feature, using 6 level set (one for each block's side), below there is a snippet of the code.

The function ls_plane() computes the distance between a point and the plane passing through a block's side, the variable cell represents the block.

The result seems to be good, the computation time for solving the problem is around 0.02 seconds, but to compute the pressure integral over the block I need up to 60.00 seconds (it depends of the block's sides orientation). Most of the computation time is used during the mls_cell.adapt().

Am I doing something wrong? how could I improve the computation time?

Thank you,

Calogero B. Rizzo

INTEGRATION_V="IM_TETRAHEDRON(3)";

INTEGRATION_P="IM_TETRAHEDRON(1)";

SIMPLEX_INTEGRATION="IM_STRUCTURED_COMPOSITE(IM_TETRAHEDRON(3),1)";

void meshCutter(const cell3D &cell) {

mls_cell.clear();

mimls_cell.clear();

getfem::level_set ls_xPos(mesh), ls_xNeg(mesh);

getfem::level_set ls_yPos(mesh), ls_yNeg(mesh);

getfem::level_set ls_zPos(mesh), ls_zNeg(mesh);

mls_cell.add_level_set(ls_xPos); mls_cell.add_level_set(ls_xNeg);

mls_cell.add_level_set(ls_yPos); mls_cell.add_level_set(ls_yNeg);

mls_cell.add_level_set(ls_zPos); mls_cell.add_level_set(ls_zNeg);

getfem::pintegration_method ppi_v = getfem::int_method_descriptor(INTEGRATION_V);

mimls_cell.set_integration_method(mesh.convex_index(), ppi_v);

getfem::pintegration_method simp_ppi = getfem::int_method_descriptor(SIMPLEX_INTEGRATION);

mimls_cell.set_simplex_im(simp_ppi); 

ls_xPos.reinit();

for(size_type d = 0; d<ls_xPos.get_mesh_fem().nb_basic_dof(); ++d) {

ls_xPos.values(0)[d] = ls_plane(ls_xPos.get_mesh_fem().point_of_basic_dof(d), cell.xPos());

}

ls_xPos.touch();

ls_yPos.reinit();

for(size_type d = 0; d<ls_yPos.get_mesh_fem().nb_basic_dof(); ++d) {

ls_yPos.values(0)[d] = ls_plane(ls_yPos.get_mesh_fem().point_of_basic_dof(d), cell.yPos());

}

ls_yPos.touch();

ls_zPos.reinit();

for(size_type d = 0; d<ls_zPos.get_mesh_fem().nb_basic_dof(); ++d) {

ls_zPos.values(0)[d] = ls_plane(ls_zPos.get_mesh_fem().point_of_basic_dof(d), cell.zPos());

}

ls_zPos.touch();

ls_xNeg.reinit();

for(size_type d = 0; d<ls_xNeg.get_mesh_fem().nb_basic_dof(); ++d) {

ls_xNeg.values(0)[d] = ls_plane(ls_xNeg.get_mesh_fem().point_of_basic_dof(d), cell.xNeg());

}

ls_xNeg.touch();

ls_yNeg.reinit();

for(size_type d = 0; d<ls_yNeg.get_mesh_fem().nb_basic_dof(); ++d) {

ls_yNeg.values(0)[d] = ls_plane(ls_yNeg.get_mesh_fem().point_of_basic_dof(d), cell.yNeg());

}

ls_yNeg.touch();

ls_zNeg.reinit();

for(size_type d = 0; d<ls_zNeg.get_mesh_fem().nb_basic_dof(); ++d) {

ls_zNeg.values(0)[d] = ls_plane(ls_zNeg.get_mesh_fem().point_of_basic_dof(d), cell.zNeg());

}

ls_zNeg.touch();

mls_cell.adapt();

mimls_cell.set_level_set_boolean_operations("a*b*c*d*e*f");

mimls_cell.adapt();

}