89 {
91 for (auto &f : recvbuf) {
93 }
95 }
96
98}
99
101 const std::string block_name, int dim) {
103
106
107 std::regex((boost::format("%s(.*)") % block_name).str())
108
109 );
110
111 for (auto bc : bcs) {
114 faces, true),
115 "get meshset ents");
117 }
118
120};
121
123 const std::string block_name, int dim) {
124 std::map<std::string, Range>
r;
125
128
129 std::regex((boost::format("%s(.*)") % block_name).str())
130
131 );
132
133 for (auto bc : bcs) {
136 faces, true),
137 "get meshset ents");
138 r[bc->getName()] = faces;
139 }
140
142}
143
145 const unsigned int cubit_bc_type) {
148 CHKERR mesh_mng->getMeshset(ms_id, cubit_bc_type, meshset);
149 return meshset;
150};
151
152static auto save_range(moab::Interface &moab,
const std::string name,
153 const Range r, std::vector<Tag> tags = {}) {
156 CHKERR moab.add_entities(*out_meshset, r);
158 CHKERR moab.write_file(name.c_str(),
"VTK",
"", out_meshset->get_ptr(), 1,
159 tags.data(), tags.size());
160 } else {
161 MOFEM_LOG(
"SELF", Sev::warning) <<
"Empty range for " << name;
162 }
164};
165
168 ParallelComm *pcomm =
171 PSTATUS_SHARED | PSTATUS_MULTISHARED,
172 PSTATUS_NOT, -1, &boundary_ents),
173 "filter_pstatus");
174 return boundary_ents;
175};
176
179 ParallelComm *pcomm =
181 CHK_MOAB_THROW(pcomm->filter_pstatus(skin, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1,
182 &owner_ents),
183 "filter_pstatus");
184 return owner_ents;
185};
186
190 CHK_MOAB_THROW(skin.find_skin(0, body_ents,
false, skin_ents),
"find_skin");
191 return skin_ents;
192};
193
196 ParallelComm *pcomm =
199 Range crack_skin_without_bdy;
200 if (pcomm->rank() == 0) {
202 CHKERR moab.get_adjacencies(crack_faces, 1,
true, crack_edges,
203 moab::Interface::UNION);
204 auto crack_skin =
get_skin(m_field, crack_faces);
208 "get_entities_by_dimension");
209 auto body_skin =
get_skin(m_field, body_ents);
210 Range body_skin_edges;
211 CHK_MOAB_THROW(moab.get_adjacencies(body_skin, 1,
true, body_skin_edges,
212 moab::Interface::UNION),
213 "get_adjacencies");
214 crack_skin_without_bdy = subtract(crack_skin, body_skin_edges);
216 for (
auto &
m : front_edges_map) {
217 auto add_front = subtract(
m.second, crack_edges);
218 auto i = intersect(
m.second, crack_edges);
220 crack_skin_without_bdy.merge(add_front);
221 } else {
224 CHKERR moab.get_adjacencies(i_skin, 1,
true, adj_i_skin,
225 moab::Interface::UNION);
226 adj_i_skin = subtract(intersect(adj_i_skin,
m.second), crack_edges);
227 crack_skin_without_bdy.merge(adj_i_skin);
228 }
229 }
230 }
231 return send_type(m_field, crack_skin_without_bdy, MBEDGE);
232}
233
236
237 ParallelComm *pcomm =
239
240 MOFEM_LOG(
"EP", Sev::noisy) <<
"get_two_sides_of_crack_surface";
241
242 if (!pcomm->rank()) {
243
244 auto impl = [&](auto &saids) {
246
248
249 auto get_adj = [&](auto e, auto dim) {
252 e, dim, true, adj, moab::Interface::UNION),
253 "get adj");
254 return adj;
255 };
256
257 auto get_conn = [&](auto e) {
260 "get connectivity");
261 return conn;
262 };
263
264 constexpr bool debug =
false;
267 body_ents);
268 auto body_skin =
get_skin(m_field, body_ents);
269 auto body_skin_edges = get_adj(body_skin, 1);
270
271 auto crack_skin =
272 subtract(
get_skin(m_field, crack_faces), body_skin_edges);
273 auto crack_skin_conn = get_conn(crack_skin);
274 auto crack_skin_conn_edges = get_adj(crack_skin_conn, 1);
275 auto crack_edges = get_adj(crack_faces, 1);
276 crack_edges = subtract(crack_edges, crack_skin);
277 auto all_tets = get_adj(crack_edges, 3);
278 crack_edges = subtract(crack_edges, crack_skin_conn_edges);
279 auto crack_conn = get_conn(crack_edges);
280 all_tets.merge(get_adj(crack_conn, 3));
281
286 crack_edges);
287 }
288
289 if (crack_faces.size()) {
290 auto grow = [&](
auto r) {
291 auto crack_faces_conn = get_conn(crack_faces);
293 auto size_r = 0;
294 while (size_r !=
r.size() &&
r.size() > 0) {
296 CHKERR moab.get_connectivity(r,
v,
true);
297 v = subtract(
v, crack_faces_conn);
300 moab::Interface::UNION);
301 r = intersect(r, all_tets);
302 }
304 break;
305 }
306 }
308 };
309
310 Range all_tets_ord = all_tets;
311 while (all_tets.size()) {
312 Range faces = get_adj(unite(saids.first, saids.second), 2);
313 faces = subtract(crack_faces, faces);
314 if (faces.size()) {
316 auto fit = faces.begin();
317 for (; fit != faces.end(); ++fit) {
318 tets = intersect(get_adj(
Range(*fit, *fit), 3), all_tets);
319 if (tets.size() == 2) {
320 break;
321 }
322 }
323 if (tets.empty()) {
324 break;
325 } else {
326 saids.first.insert(tets[0]);
327 saids.first = grow(saids.first);
328 all_tets = subtract(all_tets, saids.first);
329 if (tets.size() == 2) {
330 saids.second.insert(tets[1]);
331 saids.second = grow(saids.second);
332 all_tets = subtract(all_tets, saids.second);
333 }
334 }
335 } else {
336 break;
337 }
338 }
339
340 saids.first = subtract(all_tets_ord, saids.second);
341 saids.second = subtract(all_tets_ord, saids.first);
342 }
343
345 };
346
347 std::pair<Range, Range> saids;
348 if (crack_faces.size())
350 return saids;
351 }
352
353 MOFEM_LOG(
"EP", Sev::noisy) <<
"get_two_sides_of_crack_surface <- done";
354
355 return std::pair<Range, Range>();
356}
357
359
360struct SetIntegrationAtFrontVolume {
361
363 boost::shared_ptr<Range> front_edges)
365
367 int order_col, int order_data) {
369
370 constexpr bool debug =
false;
371
372 constexpr int numNodes = 4;
373 constexpr int numEdges = 6;
374 constexpr int refinementLevels = 4;
375
376 auto &m_field = fe_raw_ptr->
mField;
377 auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
378 auto fe_handle = fe_ptr->getFEEntityHandle();
379
380 auto set_base_quadrature = [&]() {
382 int rule = 2 * order_data + 1;
386 "wrong dimension");
387 }
391 }
393 auto &gauss_pts = fe_ptr->gaussPts;
394 gauss_pts.resize(4, nb_gauss_pts, false);
395 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[1], 4,
396 &gauss_pts(0, 0), 1);
397 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[2], 4,
398 &gauss_pts(1, 0), 1);
399 cblas_dcopy(nb_gauss_pts, &
QUAD_3D_TABLE[rule]->points[3], 4,
400 &gauss_pts(2, 0), 1);
402 &gauss_pts(3, 0), 1);
403 auto &data = fe_ptr->dataOnElement[
H1];
404 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4,
405 false);
406 double *shape_ptr =
407 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
408 cblas_dcopy(4 * nb_gauss_pts,
QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
409 1);
410 } else {
413 }
415 };
416
417 CHKERR set_base_quadrature();
418
420
421 auto get_singular_nodes = [&]() {
422 int num_nodes;
424 CHKERR m_field.
get_moab().get_connectivity(fe_handle, conn, num_nodes,
425 true);
426 std::bitset<numNodes> singular_nodes;
427 for (auto nn = 0; nn != numNodes; ++nn) {
429 singular_nodes.set(nn);
430 } else {
431 singular_nodes.reset(nn);
432 }
433 }
434 return singular_nodes;
435 };
436
437 auto get_singular_edges = [&]() {
438 std::bitset<numEdges> singular_edges;
439 for (int ee = 0; ee != numEdges; ee++) {
443 singular_edges.set(ee);
444 } else {
445 singular_edges.reset(ee);
446 }
447 }
448 return singular_edges;
449 };
450
451 auto set_gauss_pts = [&](auto &ref_gauss_pts) {
453 fe_ptr->gaussPts.swap(ref_gauss_pts);
454 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
455 auto &data = fe_ptr->dataOnElement[
H1];
456 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).resize(nb_gauss_pts, 4);
457 double *shape_ptr =
458 &*data->dataOnEntities[MBVERTEX][0].getN(
NOBASE).data().begin();
460 &fe_ptr->gaussPts(1, 0), &fe_ptr->gaussPts(2, 0),
461 nb_gauss_pts);
463 };
464
465 auto singular_nodes = get_singular_nodes();
466 if (singular_nodes.count()) {
467 auto it_map_ref_coords =
mapRefCoords.find(singular_nodes.to_ulong());
469 CHKERR set_gauss_pts(it_map_ref_coords->second);
471 } else {
472
473 auto refine_quadrature = [&]() {
475
476 const int max_level = refinementLevels;
478
479 moab::Core moab_ref;
480 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
482 for (int nn = 0; nn != 4; nn++)
483 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
484 CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
487 {
488 Range tets(tet, tet);
491 tets, 1, true, edges, moab::Interface::UNION);
494 }
495
496 Range nodes_at_front;
497 for (int nn = 0; nn != numNodes; nn++) {
498 if (singular_nodes[nn]) {
500 CHKERR moab_ref.side_element(tet, 0, nn, ent);
501 nodes_at_front.insert(ent);
502 }
503 }
504
505 auto singular_edges = get_singular_edges();
506
508 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
509 for (int ee = 0; ee != numEdges; ee++) {
510 if (singular_edges[ee]) {
512 CHKERR moab_ref.side_element(tet, 1, ee, ent);
513 CHKERR moab_ref.add_entities(meshset, &ent, 1);
514 }
515 }
516
517
519 for (int ll = 0; ll != max_level; ll++) {
522 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(ll),
524 edges);
526 CHKERR moab_ref.get_adjacencies(
527 nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
528 ref_edges = intersect(ref_edges, edges);
530 CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents,
true);
531 ref_edges = intersect(ref_edges, ents);
534 ->getEntitiesByTypeAndRefLevel(
536 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
540 ->updateMeshsetByEntitiesChildren(meshset,
542 meshset, MBEDGE, true);
543 }
544
545
548 ->getEntitiesByTypeAndRefLevel(
BitRefLevel().set(max_level),
550 tets);
551
554 }
555
557 int tt = 0;
558 for (Range::iterator tit = tets.begin(); tit != tets.end();
559 tit++, tt++) {
560 int num_nodes;
562 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes,
false);
563 CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
564 }
565
566 auto &data = fe_ptr->dataOnElement[
H1];
567 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
568 MatrixDouble ref_gauss_pts(4, nb_gauss_pts * ref_coords.size1());
570 data->dataOnEntities[MBVERTEX][0].getN(
NOBASE);
571 int gg = 0;
572 for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
573 double *tet_coords = &ref_coords(tt, 0);
575 det *= 6;
576 for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
577 for (
int dd = 0;
dd != 3;
dd++) {
578 ref_gauss_pts(dd, gg) =
579 shape_n(ggg, 0) * tet_coords[3 * 0 +
dd] +
580 shape_n(ggg, 1) * tet_coords[3 * 1 +
dd] +
581 shape_n(ggg, 2) * tet_coords[3 * 2 +
dd] +
582 shape_n(ggg, 3) * tet_coords[3 * 3 +
dd];
583 }
584 ref_gauss_pts(3, gg) = fe_ptr->gaussPts(3, ggg) * det;
585 }
586 }
587
588 mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
590
static auto send_type(MoFEM::Interface &m_field, Range r, const EntityType type)
static auto get_block_meshset(MoFEM::Interface &m_field, const int ms_id, const unsigned int cubit_bc_type)
static auto get_two_sides_of_crack_surface(MoFEM::Interface &m_field, Range crack_faces)
static auto get_range_from_block_map(MoFEM::Interface &m_field, const std::string block_name, int dim)
static auto filter_owners(MoFEM::Interface &m_field, Range skin)
static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin)
static auto get_skin(MoFEM::Interface &m_field, Range body_ents)
static auto get_crack_front_edges(MoFEM::Interface &m_field, Range crack_faces)
Range get_range_from_block(MoFEM::Interface &m_field, const std::string block_name, int dim)
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MOFEM_LOG(channel, severity)
Log.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
auto ent_form_type_and_id(const EntityType type, const EntityID id)
get entity handle from type and id
#define QUAD_3D_TABLE_SIZE
static QUAD *const QUAD_3D_TABLE[]
FTensor::Index< 'm', 3 > m
static PetscBool setSingularity
MoFEMErrorCode operator()(ForcesAndSourcesCore *fe_raw_ptr, int order_row, int order_col, int order_data)
boost::shared_ptr< Range > frontNodes
static std::map< long int, MatrixDouble > mapRefCoords
SetIntegrationAtFrontVolume(boost::shared_ptr< Range > front_nodes, boost::shared_ptr< Range > front_edges)
boost::shared_ptr< Range > frontEdges
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Deprecated interface functions.
structure to get information form mofem into EntitiesFieldData
Mesh refinement interface.
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.