318 {
320
322
324
325 try {
327 } catch (const out_of_range &e) {
328 SETERRQ(
330 "Generation of reference elements for type <%s> is not implemented",
331 moab::CN::EntityTypeName(
type));
332 }
333
334 auto set_gauss_pts = [&](auto &level_gauss_pts_on_ref_mesh, auto &level_ref,
335 auto &level_shape_functions,
336
337 auto start_vert_handle, auto start_ele_handle,
338 auto &verts_array, auto &conn, auto &ver_count,
339 auto &ele_count
340
341 ) {
343
345 level = std::min(level, level_gauss_pts_on_ref_mesh.size() - 1);
346
347 auto &level_ref_gauss_pts = level_gauss_pts_on_ref_mesh[level];
348 auto &level_ref_ele = level_ref[level];
349 auto &shape_functions = level_shape_functions[level];
350 E::gaussPts.resize(level_ref_gauss_pts.size1(), level_ref_gauss_pts.size2(),
351 false);
352 noalias(E::gaussPts) = level_ref_gauss_pts;
353
354 const auto fe_ent = E::numeredEntFiniteElementPtr->getEnt();
355 auto get_fe_coords = [&]() {
357 int num_nodes;
359 E::mField.get_moab().get_connectivity(fe_ent, conn, num_nodes, true),
360 "error get connectivity");
363 E::mField.get_moab().get_coords(conn, num_nodes, &*coords.begin()),
364 "error get coordinates");
365 return coords;
366 };
367
368 auto coords = get_fe_coords();
369
370 const int num_nodes = level_ref_gauss_pts.size2();
372
375 &*shape_functions.data().begin());
377 &verts_array[0][ver_count], &verts_array[1][ver_count],
378 &verts_array[2][ver_count]);
379 for (int gg = 0; gg != num_nodes; ++gg, ++ver_count) {
380
382
383 auto set_float_precision = [](const double x) {
384 if (std::abs(x) < std::numeric_limits<float>::epsilon())
385 return 0.;
386 else
387 return x;
388 };
389
391 auto t_ele_coords = getFTensor1FromArray<3, 3>(coords);
392 for (
int nn = 0; nn != CN::VerticesPerEntity(
type); ++nn) {
393 t_coords(
i) += t_n * t_ele_coords(
i);
394 ++t_ele_coords;
395 ++t_n;
396 }
397
398 for (auto ii : {0, 1, 2})
399 t_coords(ii) = set_float_precision(t_coords(ii));
400
401 ++t_coords;
402 }
403
405 int def_in_the_loop = -1;
407 "NB_IN_THE_LOOP", 1, MB_TYPE_INTEGER,
th, MB_TAG_CREAT | MB_TAG_SPARSE,
408 &def_in_the_loop);
409
411 const int num_el = level_ref_ele.size1();
412 const int num_nodes_on_ele = level_ref_ele.size2();
413 auto start_e = start_ele_handle + ele_count;
415 for (auto tt = 0; tt != level_ref_ele.size1(); ++tt, ++ele_count) {
416 for (int nn = 0; nn != num_nodes_on_ele; ++nn) {
417 conn[num_nodes_on_ele * ele_count + nn] =
419 }
420 }
421
422 const int n_in_the_loop = E::nInTheLoop;
424 &n_in_the_loop);
426
428 };
429
431
432 ref_ele->levelGaussPtsOnRefMesh, ref_ele->levelRef,
433 ref_ele->levelShapeFunctions,
434
435 ref_ele->startingVertEleHandle, ref_ele->startingEleHandle,
436 ref_ele->verticesOnEleArrays, ref_ele->eleConn, ref_ele->countVertEle,
437 ref_ele->countEle
438
439 );
440
442};
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
FTensor::Index< 'i', SPACE_DIM > i
UBlasVector< double > VectorDouble
auto type_from_handle(const EntityHandle h)
get type from entity handle
virtual MoFEMErrorCode transferTags()