Working to demo 3d XFEM crack propagation of a steel compact specimen top-half (gmsh import).
It compiles fine but the executable throws a "level is zero" error which I'm uncertain of the source of this error. Code attached below.
/address@hidden compact.cc
@brief three-dimensional XFEM crack propagation.
*/
#include "getfem/getfem_model_solvers.h"
#include "getfem/getfem_export.h"
#include "gmm/gmm.h"
#include "getfem/getfem_import.h"
#include "getfem/getfem_assembling.h"
#include "getfem/getfem_generic_assembly.h"
#include "getfem/getfem_mesh_im_level_set.h"
#include "getfem/getfem_mesh_fem_level_set.h"
#include "getfem/getfem_mesh_fem_product.h"
#include "getfem/getfem_mesh_fem_global_function.h"
#include "getfem/getfem_mesh_fem_sum.h"
#include <string>
#include <cstring>
#include <sstream>
#include <iostream>
using std::endl; using std::cout; using std::cerr;
using std::ends; using std::cin;
using bgeot::size_type;
using bgeot::base_node;
using bgeot::base_small_vector;
typedef getfem::model_real_plain_vector plain_vector;
using bgeot::scalar_type; /* = double */
int main(void) {
double E = 21E6; // Young Modulus (N/cm^2)
double nu = 0.3; // Poisson ratio
double clambda = E*nu/((1+nu)*(1-2*nu));
double cmu = E/(2*(1+nu));
double clambdastar = 2*clambda*cmu/(clambda+2*cmu);
double residual = 1E-9;
getfem::mesh mesh;
getfem::import_mesh("gmsh:v4.gmsh.msh",mesh);
getfem::mesh_fem mf(mesh, 3);
getfem::mesh_fem mfvm(mesh, 1);
mf.set_finite_element(mesh.convex_index(), getfem::fem_descriptor("FEM_PK(3,2)"));
mfvm.set_classical_discontinuous_finite_element(1);
getfem::level_set ls(mesh, 1, true);
getfem::mesh_level_set mls(mesh);
mls.add_level_set(ls);
mls.adapt();
std::vector<getfem::pglobal_function> vfunc(6);
for (unsigned i = 0; i < vfunc.size(); ++i){
getfem::abstract_xy_function *s =
new getfem::crack_singular_xy_function(i);
getfem::abstract_xy_function *c =
new getfem::cutoff_xy_function(0, .2,.01,.49);
s = new getfem::product_of_xy_functions(*s, *c);
vfunc[i] = getfem::global_function_on_level_set(ls, *s);
}
getfem::mesh_fem_level_set mfls(mls, mf);
getfem::mesh_fem_global_function mf_sing(mesh,3);
mf_sing.set_functions(vfunc);
getfem::mesh_fem_product mf_asympt(mfls, mf_sing);
getfem::mesh_fem_product mfu(mf_asympt);
getfem::mesh_im_level_set mim(mls, getfem::mesh_im_level_set::INTEGRATE_ALL, getfem::int_method_descriptor("IM_TETRAHEDRON(6)"), getfem::int_method_descriptor("IM_QUASI_POLAR(IM_TETRAHEDRON(8), 2)"));
getfem::model md;
md.add_fem_variable("u", mfu);
md.add_initialized_scalar_data("cmu", cmu);
md.add_initialized_scalar_data("clambdastar", clambdastar);
getfem::add_isotropic_linearized_elasticity_brick
(md, mim, "u", "clambdastar", "cmu");
getfem::add_normal_Dirichlet_condition_with_multipliers(md,mim,"u", 1, 109);
getfem::add_source_term_brick(md,mim,"u","[0,1640,0]", 108);
gmm::iteration iter(residual, 1, 50);
getfem::standard_solve(md, iter, getfem::rselect_linear_solver(md,"superlu"));
getfem::standard_solve(md, iter);
plain_vector U(mfu.nb_dof()); gmm::copy(md.real_variable("u"), U);
plain_vector VM(mfvm.nb_dof());
getfem::compute_isotropic_linearized_Von_Mises_or_Tresca
(md, "u", "clamabdastar", "cmu", mfvm, VM, false);
getfem::vtk_export exp2("d_w_vm.vtk", false);
exp2.exporting(mfu);
exp2.write_point_data(mfu, U, "displacement");
exp2.write_point_data(mfvm, VM, "Von Mises stress");
}
Thank you for any help in advance,
Andy