v0.13.1
CrackPropagation.cpp
Go to the documentation of this file.
1/** \file CrackPropagation.cpp
2 */
3
4/* This file is part of MoFEM.
5 * MoFEM is free software: you can redistribute it and/or modify it under
6 * the terms of the GNU Lesser General Public License as published by the
7 * Free Software Foundation, either version 3 of the License, or (at your
8 * option) any later version.
9 *
10 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
13 * License for more details.
14 *
15 * You should have received a copy of the GNU Lesser General Public
16 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
17
18#ifdef WITH_TETGEN
19
20#include <tetgen.h>
21#ifdef REAL
22#undef REAL
23#endif
24
25#endif
26
27#include <moab/SpatialLocator.hpp>
28#include <moab/ElemEvaluator.hpp>
29
30#include <MoFEM.hpp>
31using namespace MoFEM;
32
33#include <cholesky.hpp>
34
35#include <BasicFiniteElements.hpp>
36#include <Mortar.hpp>
37#include <Hooke.hpp>
38#include <AnalyticalFun.hpp>
39
40// TODO implementation in ComplexConstArea is obsolete, not follows naming
41// convention and code is unnecessarily complicated and in places difficult to
42// follow.
43#include <ComplexConstArea.hpp>
44#include <ConstantArea.hpp>
45
46extern "C" {
47#include <phg-quadrule/quad.h>
48}
49
50#include <NeoHookean.hpp>
51#include <Hooke.hpp>
52#include <Smoother.hpp>
55#include <Mortar.hpp>
56
57#include <CrackFrontElement.hpp>
58#include <MWLS.hpp>
60#include <CrackPropagation.hpp>
61#include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
62
63#include <CPSolvers.hpp>
64#include <CPMeshCut.hpp>
65
66#include <thread>
67
68static MoFEMErrorCode snes_monitor_fields(SNES snes, PetscInt its,
69 PetscReal fgnorm,
70 PetscViewerAndFormat *vf) {
71 PetscViewer viewer = vf->viewer;
72 Vec res;
73 DM dm;
74 PetscSection s;
75 const PetscScalar *r;
76 PetscReal *lnorms, *norms;
77 PetscInt numFields, f, pStart, pEnd, p;
78
80 LogManager::setLog("PETSC");
81 MOFEM_LOG_TAG("PETSC", "SnesProp");
82
83 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 4);
84 CHKERR SNESGetFunction(snes, &res, 0, 0);
85 CHKERR SNESGetDM(snes, &dm);
86 CHKERR DMGetDefaultSection(dm, &s);
87 CHKERR PetscSectionGetNumFields(s, &numFields);
88 CHKERR PetscSectionGetChart(s, &pStart, &pEnd);
89 CHKERR PetscCalloc2(numFields, &lnorms, numFields, &norms);
90 CHKERR VecGetArrayRead(res, &r);
91 for (p = pStart; p < pEnd; ++p) {
92 for (f = 0; f < numFields; ++f) {
93 PetscInt fdof, foff, d;
94
95 CHKERR PetscSectionGetFieldDof(s, p, f, &fdof);
96 CHKERR PetscSectionGetFieldOffset(s, p, f, &foff);
97 for (d = 0; d < fdof; ++d)
98 lnorms[f] += PetscRealPart(PetscSqr(r[foff + d]));
99 }
100 }
101 CHKERR VecRestoreArrayRead(res, &r);
102 CHKERR MPIU_Allreduce(lnorms, norms, numFields, MPIU_REAL, MPIU_SUM,
103 PetscObjectComm((PetscObject)dm));
104 CHKERR PetscViewerPushFormat(viewer, vf->format);
105 CHKERR PetscViewerASCIIAddTab(viewer, ((PetscObject)snes)->tablevel);
106 CHKERR PetscViewerASCIIPrintf(viewer, "%3D SNES Function norm %14.12e\n", its,
107 (double)fgnorm);
108 for (f = 0; f < numFields; ++f) {
109 const char *field_name;
110 CHKERR PetscSectionGetFieldName(s, f, &field_name);
111 CHKERR PetscViewerASCIIPrintf(viewer, "\t%30.30s\t%14.12e\n", field_name,
112 (double)PetscSqrtReal(norms[f]));
113 }
114 CHKERR PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel);
115 CHKERR PetscViewerPopFormat(viewer);
116 CHKERR PetscFree2(lnorms, norms);
117
118 LogManager::setLog("PETSC");
120}
121
122namespace po = boost::program_options;
123
124namespace FractureMechanics {
125
127 boost::shared_ptr<WrapMPIComm> moab_comm_wrap) {
129 std::vector<ParallelComm *> list_pcomms;
130 ParallelComm::get_all_pcomm(&moab, list_pcomms);
131 if (list_pcomms[MYPCOMM_INDEX]) {
132 CHKERR moab.tag_delete(list_pcomms[MYPCOMM_INDEX]->pstatus_tag());
133 CHKERR moab.tag_delete(list_pcomms[MYPCOMM_INDEX]->sharedps_tag());
134 CHKERR moab.tag_delete(list_pcomms[MYPCOMM_INDEX]->sharedhs_tag());
135 CHKERR moab.tag_delete(list_pcomms[MYPCOMM_INDEX]->sharedp_tag());
136 CHKERR moab.tag_delete(list_pcomms[MYPCOMM_INDEX]->sharedh_tag());
137 }
138 for (auto p : list_pcomms) {
139 delete p;
140 }
141 int pcomm_id = MYPCOMM_INDEX;
142 list_pcomms[MYPCOMM_INDEX] =
143 new ParallelComm(&moab, moab_comm_wrap->get_comm(), &pcomm_id);
145}
146
148 moab::Interface &moab_tmp,
149 const int from_proc, Range &entities,
150 const bool adjacencies, const bool tags) {
151 ErrorCode result = MB_SUCCESS;
152 int success;
153 int buff_size;
154
155 auto add_verts = [&](Range &sent_ents) {
156 // Get the verts adj to these entities, since we'll have to send those
157 // too
158
159 // First check sets
160 std::pair<Range::const_iterator, Range::const_iterator> set_range =
161 sent_ents.equal_range(MBENTITYSET);
162 ErrorCode result = MB_SUCCESS, tmp_result;
163 for (Range::const_iterator rit = set_range.first; rit != set_range.second;
164 ++rit) {
165 tmp_result = moab_tmp.get_entities_by_type(*rit, MBVERTEX, sent_ents);
166 CHK_MOAB_THROW(tmp_result, "Failed to get contained verts");
167 }
168
169 // Now non-sets
170 Range tmp_ents;
171 std::copy(sent_ents.begin(), set_range.first, range_inserter(tmp_ents));
172 result = moab_tmp.get_adjacencies(tmp_ents, 0, false, sent_ents,
173 moab::Interface::UNION);
174 CHK_MOAB_THROW(result, "Failed to get vertices adj to ghosted ents");
175
176 return result;
177 };
178
179 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab_tmp, MYPCOMM_INDEX);
180 auto &procConfig = pcomm->proc_config();
181 const int MAX_BCAST_SIZE = (1 << 28);
182
183 ParallelComm::Buffer buff(ParallelComm::INITIAL_BUFF_SIZE);
184 buff.reset_ptr(sizeof(int));
185 if ((int)procConfig.proc_rank() == from_proc) {
186
187 result = add_verts(entities);
188 CHK_MOAB_THROW(result, "Failed to add adj vertices");
189
190 buff.reset_ptr(sizeof(int));
191 result = pcomm->pack_buffer(entities, adjacencies, tags, false, -1, &buff);
192 CHK_MOAB_THROW(result,
193 "Failed to compute buffer size in broadcast_entities");
194 buff.set_stored_size();
195 buff_size = buff.buff_ptr - buff.mem_ptr;
196 }
197
198 success =
199 MPI_Bcast(&buff_size, 1, MPI_INT, from_proc, procConfig.proc_comm());
200 if (MPI_SUCCESS != success) {
201 THROW_MESSAGE("MPI_Bcast of buffer size failed");
202 }
203
204 if (!buff_size) // No data
205 return MB_SUCCESS;
206
207 if ((int)procConfig.proc_rank() != from_proc)
208 buff.reserve(buff_size);
209
210 size_t offset = 0;
211 while (buff_size) {
212 int sz = std::min(buff_size, MAX_BCAST_SIZE);
213 success = MPI_Bcast(buff.mem_ptr + offset, sz, MPI_UNSIGNED_CHAR, from_proc,
214 procConfig.proc_comm());
215 if (MPI_SUCCESS != success) {
216 THROW_MESSAGE("MPI_Bcast of buffer failed");
217 }
218
219 offset += sz;
220 buff_size -= sz;
221 }
222
223 std::vector<std::vector<EntityHandle>> dum1a, dum1b;
224 std::vector<std::vector<int>> dum1p;
225 std::vector<EntityHandle> dum2, dum4;
226 std::vector<unsigned int> dum3;
227 buff.reset_ptr(sizeof(int));
228 if ((int)procConfig.proc_rank() == from_proc) {
229 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
230 result = pcomm->unpack_buffer(buff.buff_ptr, false, from_proc, -1, dum1a,
231 dum1b, dum1p, dum2, dum2, dum3, dum4);
232 } else {
233 result = pcomm->unpack_buffer(buff.buff_ptr, false, from_proc, -1, dum1a,
234 dum1b, dum1p, dum2, dum2, dum3, dum4);
235 }
236 CHK_MOAB_THROW(result, "Failed to unpack buffer in broadcast_entities");
237 std::copy(dum4.begin(), dum4.end(), range_inserter(entities));
238
239 return MB_SUCCESS;
240};
241
246 ierr = getOptions();
247 CHKERRABORT(PETSC_COMM_WORLD, ierr);
248 }
249
252 ierr =
253 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Fixing nodes options", "none");
254 CHKERRG(ierr);
255 CHKERR PetscOptionsScalar("-gc_fix_threshold",
256 "fix nodes below given threshold", "",
257 gcFixThreshold, &gcFixThreshold, PETSC_NULL);
258 CHKERR PetscPrintf(PETSC_COMM_WORLD,
259 "### Input parameter: -gc_fix_threshold %6.4e \n",
261
262 ierr = PetscOptionsEnd();
263 CHKERRG(ierr);
265 }
266
267 MoFEMErrorCode operator()(Range &fix_material_nodes,
268 const bool fix_small_g = true,
269 const bool debug = false) const {
270 MoFEM::Interface &m_field = cP.mField;
272
273 Skinner skin(&m_field.get_moab());
274 Range tets_level, prisms_level;
275 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
276 cP.mapBitLevel["material_domain"], BitRefLevel().set(), MBTET,
277 tets_level);
278 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
279 cP.mapBitLevel["material_domain"], BitRefLevel().set(), MBPRISM,
280 prisms_level);
281
282 prisms_level = intersect(prisms_level, cP.contactElements);
283 tets_level.merge(prisms_level);
284
285 Range mesh_skin;
286 CHKERR skin.find_skin(0, tets_level, false, mesh_skin);
287
288 mesh_skin = mesh_skin.subset_by_type(MBTRI);
289
290 Range faces_to_remove, contact_faces;
291 faces_to_remove = unite(cP.crackFaces, cP.bodySkin);
292 contact_faces = unite(cP.contactMasterFaces, cP.contactSlaveFaces);
293 faces_to_remove = subtract(faces_to_remove, contact_faces);
294
295 mesh_skin = subtract(mesh_skin, faces_to_remove);
296
297 CHKERR m_field.get_moab().get_connectivity(mesh_skin, fix_material_nodes);
298 // fix fixed edges nodes
299 Range fixed_edges = cP.getInterface<CPMeshCut>()->getFixedEdges();
300 Range fixed_edges_nodes;
301 CHKERR m_field.get_moab().get_connectivity(fixed_edges, fixed_edges_nodes,
302 true);
303
304 Range crack_faces_and_fixed_edges_nodes;
305 CHKERR m_field.get_moab().get_connectivity(
306 cP.crackFaces, crack_faces_and_fixed_edges_nodes, true);
307 crack_faces_and_fixed_edges_nodes =
308 intersect(crack_faces_and_fixed_edges_nodes, fixed_edges_nodes);
309
310 Range crack_faces_edges;
311 CHKERR m_field.get_moab().get_adjacencies(
312 cP.crackFaces, 1, false, crack_faces_edges, moab::Interface::UNION);
313 Range not_fixed_edges = intersect(crack_faces_edges, fixed_edges);
314 Range not_fixed_edges_nodes;
315 CHKERR m_field.get_moab().get_connectivity(not_fixed_edges,
316 not_fixed_edges_nodes, true);
317 crack_faces_and_fixed_edges_nodes =
318 subtract(crack_faces_and_fixed_edges_nodes, not_fixed_edges_nodes);
319 fix_material_nodes.merge(crack_faces_and_fixed_edges_nodes);
320
321 Range not_fixed_edges_skin;
322 CHKERR skin.find_skin(0, not_fixed_edges, false, not_fixed_edges_skin);
323 fix_material_nodes.merge(
324 subtract(not_fixed_edges_skin, cP.crackFrontNodes));
325 if (debug && !m_field.get_comm_rank()) {
326 EntityHandle meshset;
327 CHKERR m_field.get_moab().create_meshset(MESHSET_SET, meshset);
328 CHKERR m_field.get_moab().add_entities(meshset, mesh_skin);
329 CHKERR m_field.get_moab().write_file("fixed_faces.vtk", "VTK", "",
330 &meshset, 1);
331 CHKERR m_field.get_moab().delete_entities(&meshset, 1);
332 }
333
334 // fix corner nodes
335 fix_material_nodes.merge(cP.getInterface<CPMeshCut>()->getCornerNodes());
336
337 // fix nodes on kinematic and static boundary conditions
338 auto bc_fix_nodes = [&m_field,
339 &fix_material_nodes](const EntityHandle meshset) {
341 Range ents;
342 CHKERR m_field.get_moab().get_entities_by_handle(meshset, ents, true);
343 Range nodes;
344 CHKERR m_field.get_moab().get_connectivity(ents, nodes);
345 fix_material_nodes.merge(nodes);
347 };
349 m_field, NODESET | DISPLACEMENTSET, it)) {
350 CHKERR bc_fix_nodes(it->getMeshset());
351 }
353 it)) {
354 CHKERR bc_fix_nodes(it->getMeshset());
355 }
356 if (!cP.isPressureAle) {
357 // fix nodes on surfaces under pressure
359 m_field, SIDESET | PRESSURESET, it)) {
360 CHKERR bc_fix_nodes(it->getMeshset());
361 }
362 }
363
364 // fix nodes on contact surfaces
365 if (cP.fixContactNodes) {
366 Range contact_nodes;
367 CHKERR m_field.get_moab().get_connectivity(cP.contactElements,
368 contact_nodes);
369 fix_material_nodes.merge(contact_nodes);
370 }
371
372 if (!cP.areSpringsAle) {
373 // fix nodes on spring nodes
375 if (bit->getName().compare(0, 9, "SPRING_BC") == 0) {
376 CHKERR bc_fix_nodes(bit->meshset);
377 }
378 }
379 }
380
381 {
382 Range contact_nodes;
383 CHKERR m_field.get_moab().get_connectivity(cP.mortarContactElements,
384 contact_nodes);
385 fix_material_nodes.merge(contact_nodes);
386 }
387
388 // fix nodes on the constrained interface
389 if (!cP.constrainedInterface.empty()) {
390 Range interface_nodes;
391 CHKERR m_field.get_moab().get_connectivity(cP.constrainedInterface,
392 interface_nodes);
393 fix_material_nodes.merge(interface_nodes);
394 }
395
396 auto add_nodes_on_opposite_side_of_prism = [&m_field, &fix_material_nodes,
397 this](bool only_contact_prisms =
398 false) {
400 Range adj_prisms;
401 CHKERR m_field.get_moab().get_adjacencies(
402 fix_material_nodes, 3, false, adj_prisms, moab::Interface::UNION);
403 adj_prisms = adj_prisms.subset_by_type(MBPRISM);
404 if (only_contact_prisms) {
405 adj_prisms = intersect(adj_prisms, cP.contactElements);
406 }
407 for (auto p : adj_prisms) {
408 int num_nodes;
409 const EntityHandle *conn;
410 CHKERR m_field.get_moab().get_connectivity(p, conn, num_nodes, false);
411 Range add_nodes;
412 for (int n = 0; n != 3; ++n) {
413 if (fix_material_nodes.find(conn[n]) != fix_material_nodes.end()) {
414 add_nodes.insert(conn[3 + n]);
415 }
416 if (fix_material_nodes.find(conn[3 + n]) !=
417 fix_material_nodes.end()) {
418 add_nodes.insert(conn[n]);
419 }
420 }
421 fix_material_nodes.merge(add_nodes);
422 }
424 };
425
426 CHKERR add_nodes_on_opposite_side_of_prism();
427
429 if (fix_small_g) {
430
431 if (!cP.mField.get_comm_rank()) {
432 auto copy_map = [&fix_material_nodes](const auto &map) {
433 std::map<EntityHandle, double> map_copy;
434 for (auto &m : map) {
435 if (fix_material_nodes.find(m.first) == fix_material_nodes.end()) {
436 map_copy[m.first] = m.second;
437 }
438 }
439 return map_copy;
440 };
441
442 const std::map<EntityHandle, double> map_j = copy_map(cP.mapJ);
443 const std::map<EntityHandle, double> map_g1 = copy_map(cP.mapG1);
444
445 auto find_max = [](const auto &map) {
446 double max = 0;
447 for (auto &m : map)
448 max = fmax(m.second, max);
449 return max;
450 };
451
452 const double max_j = find_max(map_j);
453 const double max_g1 = find_max(map_g1);
454
455 for (auto &m : map_j) {
456 if ((m.second < gcFixThreshold * max_j) &&
457 (map_g1.at(m.first) < gcFixThreshold * max_g1)) {
458 fix_material_nodes.insert(m.first);
459 } else {
461 }
462 }
463
464 if (map_j.empty() || map_g1.empty() || cP.fractionOfFixedNodes == 0) {
465 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
466 "No points to move on crack front");
467 }
468
469 cP.fractionOfFixedNodes /= map_j.size();
470 }
471
472 Vec vec;
473 CHKERR VecCreateMPI(cP.mField.get_comm(), PETSC_DETERMINE, 1, &vec);
474 if (!cP.mField.get_comm_rank()) {
475 CHKERR VecSetValue(vec, cP.mField.get_comm_rank(),
476 cP.fractionOfFixedNodes, INSERT_VALUES);
477 }
478 CHKERR VecMax(vec, PETSC_NULL, &cP.fractionOfFixedNodes);
479 CHKERR VecDestroy(&vec);
480 }
481
482 CHKERR add_nodes_on_opposite_side_of_prism(true);
483
485 }
486};
487
488const char *materials_list[] = {"HOOKE", "KIRCHHOFF", "NEOHOOKEAN",
489 "BONEHOOKE"};
490
492CrackPropagation::query_interface(boost::typeindex::type_index type_index,
493 UnknownInterface **iface) const {
495 *iface = NULL;
496
497 if (type_index == boost::typeindex::type_id<CrackPropagation>()) {
498 *iface = const_cast<CrackPropagation *>(this);
499 return 0;
500 }
501
502 if (type_index == boost::typeindex::type_id<CPSolvers>()) {
503 CHKERR cpSolversPtr->query_interface(type_index, iface);
504 return 0;
505 }
506
507 if (type_index == boost::typeindex::type_id<CPMeshCut>()) {
508 CHKERR cpMeshCutPtr->query_interface(type_index, iface);
509 return 0;
510 }
511 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
513}
514
516 const int approx_order,
517 const int geometry_order)
518 : mField(m_field), contactPostProcMoab(contactPostProcCore),
519 setSingularCoordinates(false), approxOrder(approx_order),
520 geometryOrder(geometry_order), gC(0), betaGc(0),
521 isGcHeterogeneous(PETSC_FALSE), isPartitioned(PETSC_FALSE),
522 propagateCrack(PETSC_FALSE), refAtCrackTip(0), refOrderAtTip(0),
523 residualStressBlock(-1), densityMapBlock(-1),
524 addAnalyticalInternalStressOperators(PETSC_FALSE), postProcLevel(0),
525 startStep(0), crackFrontLength(-1), crackSurfaceArea(-1), rHo0(1),
526 nBone(0), cpSolversPtr(new CPSolvers(*this)),
527 cpMeshCutPtr(new CPMeshCut(*this)) {
528
529 if (!LogManager::checkIfChannelExist("CPWorld")) {
530 auto core_log = logging::core::get();
531
532 core_log->add_sink(
533 LogManager::createSink(LogManager::getStrmWorld(), "CPWorld"));
534 core_log->add_sink(
535 LogManager::createSink(LogManager::getStrmSync(), "CPSync"));
536 core_log->add_sink(
537 LogManager::createSink(LogManager::getStrmSelf(), "CPSelf"));
538
539 LogManager::setLog("CPWorld");
540 LogManager::setLog("CPSync");
541 LogManager::setLog("CPSelf");
542
543 MOFEM_LOG_TAG("CPWorld", "CP");
544 MOFEM_LOG_TAG("CPSync", "CP");
545 MOFEM_LOG_TAG("CPSelf", "CP");
546 }
547
548 MOFEM_LOG("CPWorld", Sev::noisy) << "CPSolve created";
549
550#ifdef GIT_UM_SHA1_NAME
551 MOFEM_LOG_C("CPWorld", Sev::inform, "User module git commit id %s",
552 GIT_UM_SHA1_NAME);
553#endif
554
555 Version version;
556 getInterfaceVersion(version);
557 MOFEM_LOG("CPWorld", Sev::inform)
558 << "Fracture module version " << version.strVersion();
559
560#ifdef GIT_FM_SHA1_NAME
561 MOFEM_LOG_C("CPWorld", Sev::inform, "Fracture module git commit id %s",
562 GIT_FM_SHA1_NAME);
563#endif
564
565 ierr = registerInterface<CrackPropagation>();
566 CHKERRABORT(PETSC_COMM_SELF, ierr);
567 ierr = registerInterface<CPSolvers>();
568 CHKERRABORT(PETSC_COMM_SELF, ierr);
569 ierr = registerInterface<CPMeshCut>();
570 CHKERRABORT(PETSC_COMM_SELF, ierr);
571
572 mapBitLevel["mesh_cut"] = BitRefLevel().set(0);
573 mapBitLevel["spatial_domain"] = BitRefLevel().set(1);
574 mapBitLevel["material_domain"] = BitRefLevel().set(2);
575
576 if (!moabCommWorld)
577 moabCommWorld = boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
578}
579
581
584 ierr = PetscOptionsBegin(mField.get_comm(), "", "Fracture options", "none");
585 CHKERRQ(ierr);
586 {
587 CHKERR PetscOptionsBool("-my_propagate_crack",
588 "true if crack is propagated", "", propagateCrack,
589 &propagateCrack, NULL);
590 CHKERR PetscOptionsInt("-my_order", "approximation order", "", approxOrder,
591 &approxOrder, PETSC_NULL);
592 CHKERR PetscOptionsInt("-my_geom_order", "approximation geometry order", "",
593 geometryOrder, &geometryOrder, PETSC_NULL);
594 CHKERR PetscOptionsScalar("-my_gc", "release energy", "", gC, &gC,
595 PETSC_NULL);
596 CHKERR PetscOptionsScalar("-beta_gc", "heterogeneous gc beta coefficient",
598 CHKERR PetscOptionsBool("-my_is_partitioned", "true if mesh is partitioned",
599 "", isPartitioned, &isPartitioned, NULL);
600 CHKERR PetscOptionsInt("-my_ref", "crack tip mesh refinement level", "",
601 refAtCrackTip, &refAtCrackTip, PETSC_NULL);
602 CHKERR PetscOptionsInt("-my_ref_order",
603 "crack tip refinement approximation order level", "",
604 refOrderAtTip, &refOrderAtTip, PETSC_NULL);
605
606 MOFEM_LOG_C("CPWorld", Sev::inform, "### Input parameter: -my_order %d",
608 MOFEM_LOG_C("CPWorld", Sev::inform,
609 "### Input parameter: -my_geom_order %d", geometryOrder);
610 MOFEM_LOG_C("CPWorld", Sev::inform, "### Input parameter: -my_gc %6.4e",
611 gC);
612 MOFEM_LOG_C("CPWorld", Sev::inform,
613 "### Input parameter: -my_propagate_crack %d", propagateCrack);
614 MOFEM_LOG_C("CPWorld", Sev::inform, "### Input parameter: -my_ref %d",
616 MOFEM_LOG_C("CPWorld", Sev::inform, "### Input parameter: -my_ref_order %d",
618
619 onlyHookeFromOptions = PETSC_TRUE;
620 CHKERR PetscOptionsBool(
621 "-my_hook_elastic",
622 "if false force use of nonlinear element for hooke material", "",
624
625 isPressureAle = PETSC_TRUE;
626 CHKERR PetscOptionsBool("-my_add_pressure_ale",
627 "if set surface pressure is considered in ALE", "",
629
630 areSpringsAle = PETSC_TRUE;
631 CHKERR PetscOptionsBool("-my_add_springs_ale",
632 "if set surface springs is considered in ALE", "",
634
635 addSingularity = PETSC_TRUE;
636 CHKERR PetscOptionsBool("-my_add_singularity",
637 "if set singularity added to crack front", "",
639 MOFEM_LOG_C("CPWorld", Sev::inform,
640 "### Input parameter: -my_add_singularity %d", addSingularity);
641 {
642 PetscBool flg;
643 char file_name[255] = "mwls.med";
644 CHKERR PetscOptionsString(
645 "-my_mwls_approx_file",
646 "file with data from med file (code-aster) with radiation data", "",
647 file_name, file_name, 255, &flg);
648 if (flg == PETSC_TRUE) {
649 mwlsApproxFile = std::string(file_name);
650 MOFEM_LOG_C("CPWorld", Sev::inform,
651 "### Input parameter: -my_mwls_approx_file %s", file_name);
652 }
653 char tag_name[255] = "SIGP";
654 CHKERR PetscOptionsString("-my_internal_stress_name",
655 "name of internal stress tag", "", tag_name,
656 tag_name, 255, &flg);
657 MOFEM_LOG_C("CPWorld", Sev::inform,
658 "### Input parameter: -my_internal_stress_name %s", tag_name);
659 mwlsStressTagName = "MED_" + std::string(tag_name);
660 CHKERR PetscOptionsString("-my_eigen_stress_name",
661 "name of eigen stress tag", "", tag_name,
662 tag_name, 255, &flg);
663 MOFEM_LOG_C("CPWorld", Sev::inform,
664 "### Input parameter: -my_eigen_stress_name %s", tag_name);
665 mwlsEigenStressTagName = "MED_" + std::string(tag_name);
666 }
667 CHKERR PetscOptionsInt(
668 "-my_residual_stress_block",
669 "block in mechanical (i.e. cracked) mesh to which residual stress "
670 "and density mapping are applied",
671 "", residualStressBlock, &residualStressBlock, PETSC_NULL);
672 MOFEM_LOG_C("CPWorld", Sev::inform,
673 "### Input parameter: -my_residual_stress_block %d",
675
676 CHKERR PetscOptionsInt("-my_density_block",
677 "block to which density is mapped", "",
678 densityMapBlock, &densityMapBlock, PETSC_NULL);
679 MOFEM_LOG_C("CPWorld", Sev::inform,
680 "### Input parameter: -my_density_block %d", densityMapBlock);
681
682 char density_tag_name[255] = "RHO";
683 CHKERR PetscOptionsString("-my_density_tag_name",
684 "name of density tag (RHO)", "", density_tag_name,
685 density_tag_name, 255, PETSC_NULL);
686 MOFEM_LOG_C("CPWorld", Sev::inform,
687 "### Input parameter: -my_density_tag_name %s",
688 density_tag_name);
689 mwlsRhoTagName = std::string(density_tag_name);
690
691 {
692 PetscBool flg;
693 char file_name[255] = "add_cubit_meshsets.in";
694 CHKERR PetscOptionsString("-meshsets_config",
695 "meshsets config file name", "", file_name,
696 file_name, 255, &flg);
697 if (flg == PETSC_TRUE) {
698 ifstream f(file_name);
699 if (!f.good()) {
701 "File configuring meshsets ( %s ) cannot be open\n",
702 file_name);
703 }
704 configFile = file_name;
705 }
706 }
707 {
708 PetscBool flg;
710 CHKERR PetscOptionsEList("-material", "Material type", "", materials_list,
712 &flg);
713 MOFEM_LOG_C("CPWorld", Sev::inform, "### Input parameter: -material %s",
715 }
716 // Set paramaters with bone density like material
717 {
718 CHKERR PetscOptionsScalar("-my_rho0", "release energy", "", rHo0, &rHo0,
719 PETSC_NULL);
720 MOFEM_LOG_C("CPWorld", Sev::inform, "### Input parameter: -my_rho0 %6.4e",
721 rHo0);
722
723 CHKERR PetscOptionsScalar("-my_n_bone", "release energy", "", nBone,
724 &nBone, PETSC_NULL);
725 MOFEM_LOG_C("CPWorld", Sev::inform,
726 "### Input parameter: -my_n_bone %6.4e", nBone);
727 }
728 {
729 nbLoadSteps = 1;
730 CHKERR PetscOptionsInt("-nb_load_steps", "number of load steps", "",
731 nbLoadSteps, &nbLoadSteps, PETSC_NULL);
732 nbCutSteps = 0;
733 CHKERR PetscOptionsInt("-nb_cut_steps", "number of cut mesh steps", "",
734 nbCutSteps, &nbCutSteps, PETSC_NULL);
735 loadScale = 1;
736 CHKERR PetscOptionsScalar("-load_scale", "load scale", "", loadScale,
737 &loadScale, PETSC_NULL);
738 MOFEM_LOG_C("CPWorld", Sev::inform,
739 "### Input parameter: -load_scale %6.4e", loadScale);
740 }
741 {
742 otherSideConstrains = PETSC_FALSE;
743 CHKERR PetscOptionsBool(
744 "-other_side_constrain",
745 "Set surface constrain on other side of crack surface", "",
747 }
748 // Parameters for smoothing
749 {
750 smootherAlpha = 1;
751 smootherGamma = 0;
752 CHKERR PetscOptionsReal("-smoother_alpha", "Control mesh smoother", "",
753 smootherAlpha, &smootherAlpha, PETSC_NULL);
754 CHKERR PetscOptionsReal("-my_gamma", "Controls mesh smoother", "",
755 smootherGamma, &smootherGamma, PETSC_NULL);
756
757 MOFEM_LOG_C("CPWorld", Sev::inform,
758 "### Input parameter: -smoother_alpha %6.4e", smootherAlpha);
759 }
760 // Get Parameters for arc length
761 {
762 arcAlpha = 1;
763 arcBeta = 0;
764 arcS = 0;
765 CHKERR PetscOptionsReal("-arc_alpha", "Arc length alpha", "", arcAlpha,
766 &arcAlpha, PETSC_NULL);
767 CHKERR PetscOptionsReal("-arc_beta", "Arc length beta", "", arcBeta,
768 &arcBeta, PETSC_NULL);
769 CHKERR PetscOptionsReal("-arc_s", "Arc length step size", "", arcS, &arcS,
770 PETSC_NULL);
771 MOFEM_LOG_C("CPWorld", Sev::inform,
772 "### Input parameter: -arc_alpha %6.4e", arcAlpha);
773 MOFEM_LOG_C("CPWorld", Sev::inform,
774 "### Input parameter: -arc_beta %6.4e", arcBeta);
775 MOFEM_LOG_C("CPWorld", Sev::inform, "### Input parameter: -arc_s %6.4e",
776 arcS);
777 }
778 // Get parameters for contact
779 {
780 contactOrder = 2;
782 rValue = 1;
783 cnValue = 1;
784 ignoreContact = PETSC_FALSE;
785 fixContactNodes = PETSC_FALSE;
786 printContactState = PETSC_FALSE;
787 contactOutputIntegPts = PETSC_TRUE;
788 CHKERR PetscOptionsBool("-my_ignore_contact", "If true ignore contact",
789 "", ignoreContact, &ignoreContact, NULL);
790 CHKERR PetscOptionsBool("-my_fix_contact_nodes",
791 "If true fix contact nodes", "", fixContactNodes,
792 &fixContactNodes, NULL);
793 CHKERR PetscOptionsBool("-my_print_contact_state",
794 "If true print contact state", "",
796 CHKERR PetscOptionsBool(
797 "-my_contact_output_integ_pts",
798 "If true output data at contact integration points", "",
800 CHKERR PetscOptionsReal("-my_r_value", "Contact regularisation parameter",
801 "", rValue, &rValue, PETSC_NULL);
802 CHKERR PetscOptionsReal("-my_cn_value", "Contact augmentation parameter",
803 "", cnValue, &cnValue, PETSC_NULL);
804 CHKERR PetscOptionsInt("-my_contact_order", "contact approximation order",
805 "", contactOrder, &contactOrder, PETSC_NULL);
806 CHKERR PetscOptionsInt(
807 "-my_contact_lambda_order", "contact Lagrange multipliers order", "",
809 MOFEM_LOG_C("CPWorld", Sev::inform,
810 "### Input parameter: -my_cn_value %6.4e", cnValue);
811 MOFEM_LOG_C("CPWorld", Sev::inform,
812 "### Input parameter: -my_contact_order %d", contactOrder);
813 MOFEM_LOG_C("CPWorld", Sev::inform,
814 "### Input parameter: -my_cn_value %6.4e", cnValue);
815 MOFEM_LOG_C("CPWorld", Sev::inform,
816 "### Input parameter: -my_contact_order %d", contactOrder);
817 }
818 }
819 // Get Regine/Split/Cut options
820 {
821 doCutMesh = PETSC_FALSE;
822 CHKERR PetscOptionsBool("-cut_mesh", "If true mesh is cut by surface", "",
823 doCutMesh, &doCutMesh, NULL);
825 CHKERR PetscOptionsReal("-cut_factor", "Crack acceleration factor", "",
827 PETSC_NULL);
828 doElasticWithoutCrack = PETSC_FALSE;
829 CHKERR PetscOptionsBool(
830 "-run_elastic_without_crack", "If true mesh is cut by surface", "",
832
833 MOFEM_LOG_C("CPWorld", Sev::inform, "### Input parameter: -cut_mesh %d",
834 doCutMesh);
835 MOFEM_LOG_C("CPWorld", Sev::inform,
836 "### Input parameter: -cut_factor %6.4e",
838 MOFEM_LOG_C("CPWorld", Sev::inform,
839 "### Input parameter: -run_elastic_without_crack %d",
841 }
842 CHKERR PetscOptionsInt("-post_proc_level", "level of output files", "",
843 postProcLevel, &postProcLevel, PETSC_NULL);
844 MOFEM_LOG_C("CPWorld", Sev::verbose,
845 "### Input parameter: -post_proc_level %6.4e", postProcLevel);
846
848 CHKERR PetscOptionsBool(
849 "-add_analytical_internal_stress_operators",
850 "If true add analytical internal stress operators for mwls test", "",
853
854 solveEigenStressProblem = PETSC_FALSE;
855 CHKERR PetscOptionsBool("-solve_eigen_problem", "If true solve eigen problem",
857 NULL);
859 CHKERR PetscOptionsBool(
860 "-use_eigen_pos_simple_contact",
861 "If true use eigen positions for matching meshes contact", "",
863
864 MOFEM_LOG_C("CPWorld", Sev::inform,
865 "### Input parameter: -solve_eigen_problem %d",
867 MOFEM_LOG_C("CPWorld", Sev::inform,
868 "### Input parameter: -use_eigen_pos_simple_contact %d",
870
871 // Partitioning
873 CHKERR PetscOptionsScalar("-part_weight_power", "Partitioning weight power",
875 &partitioningWeightPower, PETSC_NULL);
876
878 CHKERR PetscOptionsBool(
879 "-calc_mwls_coeffs_every_propagation_step",
880 "If true recaulculate MWLS coefficients every propagation step", "",
882 NULL);
883
884 ierr = PetscOptionsEnd();
885 CHKERRQ(ierr);
886
888
889 CHKERRQ(ierr);
890
891 // adding mums options
892 char mumps_options[] =
893 "-mat_mumps_icntl_14 800 -mat_mumps_icntl_24 1 -mat_mumps_icntl_13 1 "
894 "-propagation_fieldsplit_0_mat_mumps_icntl_14 800 "
895 "-propagation_fieldsplit_0_mat_mumps_icntl_24 1 "
896 "-propagation_fieldsplit_0_mat_mumps_icntl_13 1 "
897 "-propagation_fieldsplit_1_mat_mumps_icntl_14 800 "
898 "-propagation_fieldsplit_1_mat_mumps_icntl_24 1 "
899 "-propagation_fieldsplit_1_mat_mumps_icntl_13 1 "
900 "-mat_mumps_icntl_20 0 "
901 "-propagation_fieldsplit_0_mat_mumps_icntl_20 0 "
902 "-propagation_fieldsplit_1_mat_mumps_icntl_20 0";
903 CHKERR PetscOptionsInsertString(NULL, mumps_options);
904
905 // Get options for interfaces
906 CHKERR getInterface<CPMeshCut>()->getOptions();
907 CHKERR getInterface<CPSolvers>()->getOptions();
908
910}
911
913 PetscBool flg;
914 enum tests {
915 ANALYTICAL_MODE_1,
916 CIRCULAR_PLATE_IN_INFINITE_BODY,
917 ANGLE_ANALYTICAL_MODE_1,
918 MWLS_INTERNAL_STRESS,
919 MWLS_INTERNAL_STRESS_ANALYTICAL,
920 CUT_CIRCLE_PLATE,
921 NOTCH_WITH_DENSITY,
922 LASTOP
923 };
924 const char *list[] = {"analytical_mode_1",
925 "circular_plate_in_infinite_body",
926 "angle_analytical_mode_1",
927 "mwls_internal_stress",
928 "mwls_internal_stress_analytical",
929 "cut_circle_plate",
930 "notch_with_density"};
931
932 int choice_value = ANALYTICAL_MODE_1;
933 double test_tol = 0;
935 CHKERR PetscOptionsBegin(mField.get_comm(), "", "Testing code", "none");
936 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-test", list, LASTOP,
937 &choice_value, &flg);
938 CHKERR PetscOptionsScalar("-test_tol", "test tolerance", "", test_tol,
939 &test_tol, PETSC_NULL);
940 ierr = PetscOptionsEnd();
941 CHKERRQ(ierr);
942
943 if (flg == PETSC_TRUE) {
944
945 if (mField.get_comm_rank() == 0) {
946
947 switch (choice_value) {
948 case ANALYTICAL_MODE_1: {
949 double ave_g1 = 0;
950 for (map<EntityHandle, double>::iterator mit = mapG1.begin();
951 mit != mapG1.end(); mit++) {
952 ave_g1 += mit->second;
953 }
954 ave_g1 /= mapG1.size();
955 // Exact solution
956 const double gc = 1e-5;
957 // Tolerance
958 const double error = fabs(ave_g1 - gc) / gc;
959 CHKERR PetscPrintf(
960 PETSC_COMM_SELF,
961 "ave_gc = %6.4e gc = %6.4e tolerance = %6.4e error = %6.4e\n",
962 ave_g1, gc, test_tol, error);
963 if (error > test_tol || error != error) {
964 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
965 "ANALYTICAL_MODE_1 test failed");
966 }
967 } break;
968 case CIRCULAR_PLATE_IN_INFINITE_BODY: {
969 const double gc = 1.1672e+01;
970 double max_error = 0;
971 for (map<EntityHandle, double>::iterator mit = mapG1.begin();
972 mit != mapG1.end(); mit++) {
973 const double g = mit->second;
974 const double error = fabs(g - gc) / gc;
975 max_error = (max_error < error) ? error : max_error;
976 CHKERR PetscPrintf(
977 PETSC_COMM_SELF,
978 "g = %6.4e gc = %6.4e tolerance = %6.4e error = %6.4e\n", g, gc,
979 test_tol, error);
980 }
981 if (max_error > test_tol || max_error != max_error)
982 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
983 "CIRCULAR_PLATE_IN_INFINITE_BODY test failed");
984
985 } break;
986 case ANGLE_ANALYTICAL_MODE_1: {
987 const double gc = 1e-5;
988 double max_error = 0;
989 for (map<EntityHandle, double>::iterator mit = mapG1.begin();
990 mit != mapG1.end(); mit++) {
991 const double g = mit->second;
992 const double error = fabs(g - gc) / gc;
993 max_error = (max_error < error) ? error : max_error;
994 CHKERR PetscPrintf(
995 PETSC_COMM_SELF,
996 "g = %6.4e gc = %6.4e tolerance = %6.4e error = %6.4e\n", g, gc,
997 test_tol, error);
998 }
999 if (max_error > test_tol || max_error != max_error) {
1000 SETERRQ2(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1001 "ANGLE_ANALYTICAL_MODE_1 test failed %6.4e > %6.4e",
1002 max_error, test_tol);
1003 }
1004 } break;
1005 case MWLS_INTERNAL_STRESS: {
1006 const double gc =
1007 1.5527e+06; // this is max g1 for big problem (not extact solution)
1008 double max_error = 0;
1009 for (map<EntityHandle, double>::iterator mit = mapG1.begin();
1010 mit != mapG1.end(); mit++) {
1011 const double g = mit->second;
1012 const double error = fabs(g - gc) / gc;
1013 max_error = (max_error < error) ? error : max_error;
1014 CHKERR PetscPrintf(
1015 PETSC_COMM_SELF,
1016 "g = %6.4e gc = %6.4e tolerance = %6.4e error = %6.4e\n", g, gc,
1017 test_tol, error);
1018 }
1019 if (max_error > test_tol || max_error != max_error) {
1020 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1021 "MWLS_INTERNAL_STRESS test failed");
1022 }
1023 } break;
1024 case MWLS_INTERNAL_STRESS_ANALYTICAL: {
1025 const double gc =
1026 1.5698e+06; // this is max g1 for big problem (not extact solution)
1027 double max_error = 0;
1028 for (map<EntityHandle, double>::iterator mit = mapG1.begin();
1029 mit != mapG1.end(); mit++) {
1030 const double g = mit->second;
1031 const double error = fabs(g - gc) / gc;
1032 max_error = (max_error < error) ? error : max_error;
1033 CHKERR PetscPrintf(
1034 PETSC_COMM_SELF,
1035 "g = %6.4e gc = %6.4e tolerance = %6.4e error = %6.4e\n", g, gc,
1036 test_tol, error);
1037 }
1038 if (max_error > test_tol || max_error != max_error) {
1039 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1040 "MWLS_INTERNAL_STRESS_ANALYTICAL test failed");
1041 }
1042 } break;
1043 case CUT_CIRCLE_PLATE: {
1044 // This is simple regression test
1045 double max_g = 0;
1046 for (auto &m : mapG1) {
1047 const double g = m.second;
1048 max_g = fmax(g, max_g);
1049 }
1050 const double gc =
1051 7.5856e-02; // this is max g1 for big problem (not extact solution)
1052 double error = fabs(max_g - gc) / gc;
1053 if (error > test_tol || error != error) {
1054 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1055 "CUT_CIRCLE_PLATE test failed");
1056 }
1057 } break;
1058 case NOTCH_WITH_DENSITY: {
1059 // This is simple regression test
1060 double max_g = 0;
1061 for (map<EntityHandle, double>::iterator mit = mapG1.begin();
1062 mit != mapG1.end(); mit++) {
1063 const double g = mit->second;
1064 max_g = fmax(g, max_g);
1065 }
1066 const double gc = 1.6642e-04; // this is max g1 for a small problem (not
1067 // extact solution)
1068 double error = fabs(max_g - gc) / gc;
1069 if (error > test_tol || error != error) {
1070 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1071 "NOTCH_WITH_DENSITY test failed");
1072 }
1073 } break;
1074 }
1075 }
1076 }
1077
1079}
1080
1082 const Range &ents) {
1084 std::ostringstream file;
1085 file << prefix << mField.get_comm_rank() << ".vtk";
1086 EntityHandle meshset_out;
1087 CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset_out);
1088 CHKERR mField.get_moab().add_entities(meshset_out, ents);
1089 CHKERR mField.get_moab().write_file(file.str().c_str(), "VTK", "",
1090 &meshset_out, 1, NULL, 0);
1091 CHKERR mField.get_moab().delete_entities(&meshset_out, 1);
1093};
1094
1096 BitRefLevel bit2, int verb,
1097 const bool debug) {
1099
1102
1103 // create map with layers with different weights for partitioning
1104 map<int, Range> layers;
1105 Range ref_nodes = crackFrontNodes;
1106 Range tets;
1107 CHKERR moab.get_adjacencies(ref_nodes, 3, false, tets,
1108 moab::Interface::UNION);
1109 CHKERR mField.getInterface<BitRefManager>()->filterEntitiesByRefLevel(
1110 bit1, BitRefLevel().set(), tets);
1111 layers[0] = tets;
1112 for (int ll = 0; ll != refOrderAtTip; ll++) {
1113
1114 CHKERR moab.get_connectivity(tets, ref_nodes, false);
1115 CHKERR moab.get_adjacencies(ref_nodes, 3, false, tets,
1116 moab::Interface::UNION);
1117 CHKERR mField.getInterface<BitRefManager>()->filterEntitiesByRefLevel(
1118 bit1, BitRefLevel().set(), tets);
1119
1120 layers[ll + 1] = subtract(tets, layers[ll]);
1121
1122 }
1123 Range tets_bit_2;
1124 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1125 bit2, BitRefLevel().set(), MBTET, tets_bit_2);
1126 int last = layers.size();
1127 auto rest_elements_in_the_bubble = subtract(tets_bit_2, layers[last - 1]);
1128 if (rest_elements_in_the_bubble.size())
1129 layers[last] = rest_elements_in_the_bubble;
1130
1131 // create tag with weights
1132 Tag th_tet_weight;
1133 rval = moab.tag_get_handle("TETS_WEIGHT", th_tet_weight);
1134 if (rval == MB_SUCCESS) {
1135 CHKERR moab.tag_delete(th_tet_weight);
1136 }
1137 int def_val = 1;
1138 CHKERR moab.tag_get_handle("TETS_WEIGHT", 1, MB_TYPE_INTEGER, th_tet_weight,
1139 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1140
1141 for (unsigned int ll = 0; ll != layers.size(); ll++) {
1142 int weight = pow(layers.size() + 2 - ll, partitioningWeightPower);
1143 CHKERR moab.tag_clear_data(th_tet_weight, layers[ll], &weight);
1144 }
1145
1146 // partition mesh
1147 Range ents;
1148 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1149 bit1, BitRefLevel().set(), MBTET, ents);
1150 if (debug) {
1151 CHKERR saveEachPart("out_ents_to_partition_", ents);
1152 }
1154 ents, 3, 0, mField.get_comm_size(), &th_tet_weight, NULL, NULL, verb,
1155 false);
1156 CHKERR moab.tag_delete(th_tet_weight);
1157
1159}
1160
1162 Range &proc_ents, const int verb,
1163 const bool debug) {
1165 Skinner skin(&mField.get_moab());
1167
1168 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
1169
1170 proc_ents.clear();
1171 Range all_proc_ents;
1172
1173 // get entities on processor
1174 Tag part_tag = pcomm->part_tag();
1175 Range tagged_sets;
1176 CHKERR mField.get_moab().get_entities_by_type_and_tag(
1177 0, MBENTITYSET, &part_tag, NULL, 1, tagged_sets, moab::Interface::UNION);
1178 for (Range::iterator mit = tagged_sets.begin(); mit != tagged_sets.end();
1179 mit++) {
1180 int part;
1181 CHKERR moab.tag_get_data(part_tag, &*mit, 1, &part);
1182 if (part == mField.get_comm_rank()) {
1183 CHKERR moab.get_entities_by_dimension(*mit, 3, proc_ents, true);
1184 CHKERR moab.get_entities_by_handle(*mit, all_proc_ents, true);
1185 }
1186 }
1187 proc_ents = intersect(proc_ents, tets);
1188
1189 // get body skin
1190 Range tets_skin;
1191 CHKERR skin.find_skin(0, tets, false, tets_skin);
1192
1193 // get skin on processor
1194 Range proc_ents_skin[4];
1195
1196 Range contact_tets;
1197 if (!contactElements.empty() && !ignoreContact && !fixContactNodes) {
1198 BitRefLevel bit2_again = mapBitLevel["material_domain"];
1199 Range prisms_level;
1200 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1201 bit2_again, BitRefLevel().set(), MBPRISM, prisms_level);
1202
1203 Range contact_prisms = intersect(prisms_level, contactElements);
1204
1205 Range contact_prisms_to_remove =
1206 intersect(prisms_level, mortarContactElements);
1207 prisms_level = subtract(prisms_level, contact_prisms_to_remove);
1208
1209 Range face_on_prism;
1210 CHKERR moab.get_adjacencies(contact_prisms, 2, false, face_on_prism,
1211 moab::Interface::UNION);
1212
1213 Range tris_on_prism = face_on_prism.subset_by_type(MBTRI);
1214
1215 Range adj_vols_to_prisms;
1216 CHKERR moab.get_adjacencies(tris_on_prism, 3, false, adj_vols_to_prisms,
1217 moab::Interface::UNION);
1218
1219 Range tet_level;
1220 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1221 bit2_again, BitRefLevel().set(), MBTET, tet_level);
1222
1223 Range contact_tets_from_vols = adj_vols_to_prisms.subset_by_type(MBTET);
1224
1225 contact_tets = intersect(tet_level, contact_tets_from_vols);
1226
1228 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
1229 Range contact_tets_container = contact_tets;
1230 // FIXME: can we delete this?
1231 // Distribute prisms elements
1232 // {
1233
1234 // Range tris_level;
1235
1236 // //Find tris of tets
1237 // CHKERR moab.get_adjacencies(contact_tets, 2, false, tris_level,
1238 // moab::Interface::UNION);
1239 // for (Range::iterator mit = tagged_sets.begin(); mit !=
1240 // tagged_sets.end();
1241 // mit++) {
1242 // Range part_tris;
1243 // CHKERR moab.get_entities_by_type(*mit, MBTRI, part_tris);
1244 // part_tris = intersect(part_tris, tris_level);
1245 // Range adj_vols;
1246 // CHKERR moab.get_adjacencies(part_tris, 3, false, adj_vols,
1247 // moab::Interface::UNION);
1248 // adj_vols = intersect(adj_vols, contact_tets_container);
1249 // contact_tets_container = subtract(contact_tets_container, adj_vols);
1250 // int part;
1251 // CHKERR moab.tag_get_data(part_tag, &*mit, 1, &part);
1252 // std::vector<int> tag(adj_vols.size(), part);
1253 // CHKERR moab.tag_set_data(part_tag, adj_vols, &*tag.begin());
1254 // CHKERR moab.add_entities(*mit, adj_vols);
1255 // }
1256 // }
1257 contactTets = contact_tets;
1258 proc_ents.merge(contact_tets);
1259
1260 Range contact_tets_vert;
1261 CHKERR moab.get_adjacencies(contact_tets, 0, false, contact_tets_vert,
1262 moab::Interface::UNION);
1263 Range contact_tets_edges;
1264 CHKERR moab.get_adjacencies(contact_tets, 1, false, contact_tets_edges,
1265 moab::Interface::UNION);
1266 Range contact_tets_faces;
1267 CHKERR moab.get_adjacencies(contact_tets, 2, false, contact_tets_faces,
1268 moab::Interface::UNION);
1269 proc_ents_skin[2].merge(contact_tets_faces);
1270 proc_ents_skin[1].merge(contact_tets_edges);
1271 proc_ents_skin[0].merge(contact_tets_vert);
1272 }
1273
1274 proc_ents_skin[3] = proc_ents;
1275 CHKERR skin.find_skin(0, proc_ents, false, proc_ents_skin[2]);
1276 // and get shared entities
1277 proc_ents_skin[2] = subtract(proc_ents_skin[2], tets_skin);
1278 // Note that all crack faces are shared
1279 proc_ents_skin[2].merge(crackFaces);
1280 proc_ents_skin[2].merge(contactSlaveFaces);
1281 proc_ents_skin[2].merge(contactMasterFaces);
1282 proc_ents_skin[2].merge(mortarContactSlaveFaces);
1283 proc_ents_skin[2].merge(mortarContactMasterFaces);
1284
1285 CHKERR moab.get_adjacencies(proc_ents_skin[2], 1, false, proc_ents_skin[1],
1286 moab::Interface::UNION);
1287 CHKERR moab.get_connectivity(proc_ents_skin[1], proc_ents_skin[0], true);
1288 proc_ents_skin[1].merge(crackFront);
1289 Range crack_faces_nodes;
1290 CHKERR moab.get_connectivity(crackFaces, crack_faces_nodes, true);
1291 proc_ents_skin[0].merge(crack_faces_nodes);
1292 // proc_ents_skin[0].merge(crackFrontNodes);
1293
1294 if (mField.check_field("LAMBDA_ARC_LENGTH")) {
1295 EntityHandle lambda_meshset = mField.get_field_meshset("LAMBDA_ARC_LENGTH");
1296 Range arc_length_vertices;
1297 CHKERR moab.get_entities_by_type(lambda_meshset, MBVERTEX,
1298 arc_length_vertices, true);
1299 if (arc_length_vertices.size() != 1) {
1301 "Should be one vertex in <LAMBDA_ARC_LENGTH> field but is %d",
1302 arc_length_vertices.size());
1303 }
1304 arcLengthVertex = arc_length_vertices[0];
1305 proc_ents_skin[0].merge(arc_length_vertices);
1306 } else {
1307 const double coords[] = {0, 0, 0};
1308 CHKERR moab.create_vertex(coords, arcLengthVertex);
1309 proc_ents_skin[0].insert(arcLengthVertex);
1310 }
1311 {
1312 Tag th_gid;
1313 CHKERR mField.get_moab().tag_get_handle(GLOBAL_ID_TAG_NAME, th_gid);
1314 int gid;
1315 ierr = mField.get_moab().get_number_entities_by_type(0, MBVERTEX, gid);
1316 CHKERR mField.get_moab().tag_set_data(th_gid, &arcLengthVertex, 1, &gid);
1317 }
1318
1319 Range all_ents_no_on_part;
1320 CHKERR moab.get_entities_by_handle(0, all_ents_no_on_part, true);
1321 all_ents_no_on_part = subtract(all_ents_no_on_part, all_proc_ents);
1322 for (int dd = 0; dd != 4; ++dd) {
1323 all_ents_no_on_part = subtract(all_ents_no_on_part, proc_ents_skin[dd]);
1324 }
1325 all_ents_no_on_part = subtract(
1326 all_ents_no_on_part, all_ents_no_on_part.subset_by_type(MBENTITYSET));
1327 CHKERR mField.remove_ents_from_finite_element(all_ents_no_on_part);
1328 CHKERR mField.remove_ents_from_field(all_ents_no_on_part);
1329
1330 {
1331 Range all_ents;
1332 CHKERR moab.get_entities_by_handle(0, all_ents);
1333 std::vector<unsigned char> pstat_tag(all_ents.size(), 0);
1334 CHKERR moab.tag_set_data(pcomm->pstatus_tag(), all_ents,
1335 &*pstat_tag.begin());
1336 }
1337
1338 CHKERR pcomm->resolve_shared_ents(0, proc_ents, 3, -1, proc_ents_skin);
1339
1340 const RefEntity_multiIndex *refined_ents_ptr;
1341 ierr = mField.get_ref_ents(&refined_ents_ptr);
1342 if (debug) {
1343 std::ostringstream file_skin;
1344 file_skin << "out_skin_" << mField.get_comm_rank() << ".vtk";
1345 EntityHandle meshset_skin;
1346 CHKERR moab.create_meshset(MESHSET_SET, meshset_skin);
1347 CHKERR moab.add_entities(meshset_skin, proc_ents_skin[2]);
1348 CHKERR moab.add_entities(meshset_skin, proc_ents_skin[1]);
1349 CHKERR moab.add_entities(meshset_skin, proc_ents_skin[0]);
1350
1351 CHKERR moab.write_file(file_skin.str().c_str(), "VTK", "", &meshset_skin,
1352 1);
1353 std::ostringstream file_owned;
1354 file_owned << "out_owned_" << mField.get_comm_rank() << ".vtk";
1355 EntityHandle meshset_owned;
1356 CHKERR moab.create_meshset(MESHSET_SET, meshset_owned);
1357 CHKERR moab.add_entities(meshset_owned, proc_ents);
1358 CHKERR moab.write_file(file_owned.str().c_str(), "VTK", "", &meshset_owned,
1359 1);
1360 CHKERR moab.delete_entities(&meshset_owned, 1);
1361 }
1362
1363 if (debug) {
1364 Range shared_ents;
1365 // Get entities shared with all other processors
1366 CHKERR pcomm->get_shared_entities(-1, shared_ents);
1367 std::ostringstream file_shared_owned;
1368 file_shared_owned << "out_shared_owned_" << mField.get_comm_rank()
1369 << ".vtk";
1370 EntityHandle meshset_shared_owned;
1371 CHKERR moab.create_meshset(MESHSET_SET, meshset_shared_owned);
1372 CHKERR moab.add_entities(meshset_shared_owned, shared_ents);
1373 CHKERR moab.write_file(file_shared_owned.str().c_str(), "VTK", "",
1374 &meshset_shared_owned, 1);
1375 CHKERR moab.delete_entities(&meshset_shared_owned, 1);
1376 }
1377
1379}
1380
1382 const int verb,
1383 const bool debug) {
1385 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
1387 bitEnts.clear();
1388 bitProcEnts.clear();
1389 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1390 bit, BitRefLevel().set(), MBTET, bitEnts);
1391
1392 if (debug) {
1393 std::ostringstream file_owned;
1394 file_owned << "out_bit_" << mField.get_comm_rank() << ".vtk";
1395 EntityHandle meshset;
1396 CHKERR moab.create_meshset(MESHSET_SET, meshset);
1397 CHKERR moab.add_entities(meshset, bitEnts.subset_by_dimension(3));
1398 CHKERR moab.write_file(file_owned.str().c_str(), "VTK", "", &meshset, 1);
1399 CHKERR moab.delete_entities(&meshset, 1);
1400 }
1401 // Distribute prisms elements
1402 {
1403 Tag part_tag = pcomm->part_tag();
1404 Range tagged_sets;
1405 CHKERR mField.get_moab().get_entities_by_type_and_tag(
1406 0, MBENTITYSET, &part_tag, NULL, 1, tagged_sets,
1407 moab::Interface::UNION);
1408 Range tris_level;
1409 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1410 bit, BitRefLevel().set(), MBTRI, tris_level);
1411 Range prisms_sides;
1412 prisms_sides.merge(oneSideCrackFaces);
1413 prisms_sides.merge(contactSlaveFaces);
1414 prisms_sides.merge(mortarContactSlaveFaces);
1415
1416 tris_level = intersect(tris_level, prisms_sides);
1417 Range prisms_level;
1418 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1419 bit, BitRefLevel().set(), MBPRISM, prisms_level);
1420
1421 for (Range::iterator mit = tagged_sets.begin(); mit != tagged_sets.end();
1422 mit++) {
1423 Range part_tris;
1424 CHKERR moab.get_entities_by_type(*mit, MBTRI, part_tris);
1425 part_tris = intersect(part_tris, tris_level);
1426 Range adj_prisms;
1427 CHKERR moab.get_adjacencies(part_tris, 3, false, adj_prisms,
1428 moab::Interface::UNION);
1429 adj_prisms = intersect(adj_prisms, prisms_level);
1430 prisms_level = subtract(prisms_level, adj_prisms);
1431
1432 int part;
1433 CHKERR moab.tag_get_data(part_tag, &*mit, 1, &part);
1434 std::vector<int> tag(adj_prisms.size(), part);
1435 CHKERR moab.tag_set_data(part_tag, adj_prisms, &*tag.begin());
1436 CHKERR moab.add_entities(*mit, adj_prisms);
1437 }
1438 }
1439
1442}
1443
1445 const BitRefLevel bit, const BitRefLevel mask, const bool proc_only,
1446 const bool build_fields, const int verb, const bool debug) {
1449
1450 auto write_file = [&](const std::string file_name, Range ents) {
1452 EntityHandle meshset;
1453 CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
1454 CHKERR mField.get_moab().add_entities(meshset, ents);
1455 if (mField.get_comm_rank() == 0) {
1456 CHKERR mField.get_moab().write_file(file_name.c_str(), "VTK", "",
1457 &meshset, 1);
1458 }
1459 CHKERR mField.get_moab().delete_entities(&meshset, 1);
1461 };
1462
1463 Range tets_level, tets_level_all;
1464 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1465 bit, mask, MBTET, tets_level_all);
1466 if (proc_only) {
1467 tets_level = intersect(bitProcEnts, tets_level_all);
1468 }
1469 // That is ok, since all crack faces are shared
1470 Range prisms_level;
1471 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1472 bit, mask, MBPRISM, prisms_level);
1473 if (proc_only) {
1474 prisms_level = intersect(prisms_level, bitProcEnts);
1475 }
1476 Range prisms_level_tris;
1477 rval = mField.get_moab().get_adjacencies(
1478 prisms_level, 2, false, prisms_level_tris, moab::Interface::UNION);
1479 // MoAB will report problems at crack front since quads are not well defined.
1480 if (rval != MB_SUCCESS && rval != MB_MULTIPLE_ENTITIES_FOUND) {
1481 CHKERRG(rval);
1482 } else {
1483 rval = MB_SUCCESS;
1484 }
1485 prisms_level_tris = intersect(prisms_level_tris, crackFaces);
1486 Range proc_ents = tets_level;
1487 proc_ents.merge(prisms_level_tris);
1488
1489 CHKERR mField.add_field("SPATIAL_POSITION", H1, AINSWORTH_LOBATTO_BASE, 3,
1490 MB_TAG_SPARSE, MF_ZERO, verb);
1492 CHKERR mField.add_field("EIGEN_SPATIAL_POSITIONS", H1,
1493 AINSWORTH_LEGENDRE_BASE, 3, MB_TAG_SPARSE, MF_ZERO,
1494 verb);
1495 }
1496 CHKERR mField.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE, 3,
1497 MB_TAG_SPARSE, MF_ZERO, verb);
1498
1499 // Set approx order
1500 Range dofs_ents;
1501 dofs_ents.merge(proc_ents);
1502 for (int dd = 1; dd != 3; ++dd) {
1503 CHKERR moab.get_adjacencies(proc_ents, dd, false, dofs_ents,
1504 moab::Interface::UNION);
1505 }
1506
1507 Range dofs_nodes;
1508 CHKERR moab.get_connectivity(proc_ents, dofs_nodes, true);
1509 dofs_nodes.merge(crackFrontNodes);
1510 dofs_ents.merge(crackFront);
1511
1512 EntityHandle spatial_meshset = mField.get_field_meshset("SPATIAL_POSITION");
1513 Range
1514 ents_to_remove_from_field; // Entities which are no part of problem any
1515 // more. Very likely partition has been chaned.
1516 CHKERR mField.get_moab().get_entities_by_handle(spatial_meshset,
1517 ents_to_remove_from_field);
1518 ents_to_remove_from_field = subtract(ents_to_remove_from_field, dofs_nodes);
1519 ents_to_remove_from_field = subtract(ents_to_remove_from_field, dofs_ents);
1520 CHKERR mField.remove_ents_from_field("SPATIAL_POSITION",
1521 ents_to_remove_from_field, verb);
1522 CHKERR mField.remove_ents_from_field("MESH_NODE_POSITIONS",
1523 ents_to_remove_from_field, verb);
1525 CHKERR mField.remove_ents_from_field("EIGEN_SPATIAL_POSITIONS",
1526 ents_to_remove_from_field, verb);
1527
1528 // Add entities to field
1529 auto add_spatial_field_ents = [&](const auto field_name) {
1532 verb);
1533 CHKERR mField.add_ents_to_field_by_type(prisms_level_tris, MBTRI,
1534 field_name, verb);
1536 field_name, verb);
1538 verb);
1539 CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
1540 field_name, verb);
1542 };
1543 CHKERR add_spatial_field_ents("SPATIAL_POSITION");
1545 CHKERR add_spatial_field_ents("EIGEN_SPATIAL_POSITIONS");
1546 }
1547
1548 Range current_field_ents;
1549 CHKERR mField.get_moab().get_entities_by_handle(spatial_meshset,
1550 current_field_ents);
1551 dofs_nodes = current_field_ents.subset_by_type(MBVERTEX);
1552 dofs_ents = subtract(current_field_ents, dofs_nodes);
1553 CHKERR mField.add_ents_to_field_by_type(dofs_nodes, MBVERTEX,
1554 "MESH_NODE_POSITIONS", verb);
1555
1556 auto set_spatial_field_order = [&](const auto field_name) {
1559 CHKERR mField.set_field_order(dofs_nodes, field_name, 1, verb);
1561 };
1562 CHKERR set_spatial_field_order("SPATIAL_POSITION");
1564 CHKERR set_spatial_field_order("EIGEN_SPATIAL_POSITIONS");
1565 }
1566
1567 CHKERR mField.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", geometryOrder,
1568 verb);
1569 CHKERR mField.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", geometryOrder,
1570 verb);
1571 CHKERR mField.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", geometryOrder,
1572 verb);
1573 CHKERR mField.set_field_order(dofs_nodes, "MESH_NODE_POSITIONS", 1, verb);
1574
1575 if (debug) {
1576 CHKERR write_file("global_order_" + std::to_string(approxOrder) + ".vtk",
1577 dofs_ents);
1578 }
1579
1580 Range non_ref_ents;
1581 non_ref_ents.merge(dofs_ents);
1582
1583 Range contact_faces;
1584 contact_faces.merge(contactSlaveFaces);
1585 contact_faces.merge(contactMasterFaces);
1586 contact_faces.merge(mortarContactMasterFaces);
1587 contact_faces.merge(mortarContactSlaveFaces);
1588
1589 // Set ho approx. at crack front
1590 if (refOrderAtTip > 0) {
1591 Range proc_ents;
1592 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1593 bit, mask, MBTET, proc_ents);
1594 map<int, Range> layers;
1595 Range ref_nodes = crackFrontNodes;
1596 for (int ll = 0; ll != refOrderAtTip; ll++) {
1597 Range tets;
1598 CHKERR moab.get_adjacencies(ref_nodes, 3, false, tets,
1599 moab::Interface::UNION);
1600 tets = intersect(tets_level_all, tets);
1601 layers[ll] = tets;
1602 Range ref_order_ents;
1603 CHKERR moab.get_adjacencies(tets, 2, false, ref_order_ents,
1604 moab::Interface::UNION);
1605 CHKERR moab.get_adjacencies(tets, 1, false, ref_order_ents,
1606 moab::Interface::UNION);
1607 layers[ll].merge(ref_order_ents);
1608 CHKERR moab.get_connectivity(tets, ref_nodes, false);
1609
1610 if (!contact_faces.empty()) {
1611 Range contact_layer = intersect(layers[ll], contact_faces);
1612 for (Range::iterator tit = contact_layer.begin();
1613 tit != contact_layer.end(); tit++) {
1614 Range contact_prisms, contact_tris;
1615 CHKERR moab.get_adjacencies(&*tit, 1, 3, false, contact_prisms,
1616 moab::Interface::UNION);
1617 contact_prisms = contact_prisms.subset_by_type(MBPRISM);
1618
1619 Range simple_contact_elements;
1620 if (!contactElements.empty()) {
1621 simple_contact_elements.merge(contactElements);
1622 simple_contact_elements =
1623 intersect(contact_prisms, simple_contact_elements);
1624 if (simple_contact_elements.size() != 1) {
1625 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1626 "Expecting 1 adjacent contact prism, but %d found",
1627 simple_contact_elements.size());
1628 }
1629 CHKERR moab.get_adjacencies(simple_contact_elements, 2, false,
1630 contact_tris, moab::Interface::UNION);
1631 contact_tris = contact_tris.subset_by_type(MBTRI);
1632 if (contact_tris.size() != 2) {
1633 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1634 "Expecting 2 adjacent contact tris, but %d found",
1635 contact_tris.size());
1636 }
1637 }
1638
1639 Range mortar_contact_elements;
1640 if (!mortarContactElements.empty() &&
1641 simple_contact_elements.empty()) {
1642 mortar_contact_elements.merge(mortarContactElements);
1643 mortar_contact_elements =
1644 intersect(contact_prisms, mortar_contact_elements);
1645 if (mortar_contact_elements.empty()) {
1646 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1647 "Expecting adjacent mortar contact prism, but %d found",
1648 mortar_contact_elements.size());
1649 }
1650 CHKERR moab.get_adjacencies(mortar_contact_elements, 2, false,
1651 contact_tris, moab::Interface::UNION);
1652 contact_tris = contact_tris.subset_by_type(MBTRI);
1653 if (contact_tris.size() < 2) {
1654 SETERRQ1(
1655 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1656 "Expecting 2 or more adjacent contact tris, but %d found",
1657 contact_tris.size());
1658 }
1659 }
1660
1661 contact_prisms.merge(simple_contact_elements);
1662 contact_prisms.merge(mortar_contact_elements);
1663
1664 if (contact_prisms.empty()) {
1665 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1666 "Expecting adjacent mortar contact prism, but %d found",
1667 contact_prisms.size());
1668 }
1669
1670 layers[ll].merge(contact_tris);
1671 Range contact_edges, contact_nodes;
1672 CHKERR moab.get_adjacencies(contact_tris, 1, false, contact_edges,
1673 moab::Interface::UNION);
1674 layers[ll].merge(contact_edges);
1675 CHKERR moab.get_connectivity(contact_tris, contact_nodes, false);
1676 ref_nodes.merge(contact_nodes);
1677 }
1678 }
1679 }
1680 int poo = 1;
1681 for (map<int, Range>::reverse_iterator mit = layers.rbegin();
1682 mit != layers.rend(); mit++, poo++) {
1683 mit->second = intersect(mit->second, dofs_ents);
1684 CHKERR mField.set_field_order(mit->second, "SPATIAL_POSITION",
1685 approxOrder + poo, verb);
1686 non_ref_ents = subtract(non_ref_ents, mit->second);
1687
1688 if (debug) {
1689 CHKERR write_file("layer_order_" + std::to_string(approxOrder + poo) +
1690 ".vtk",
1691 mit->second);
1692 }
1693 }
1694 }
1695
1696 if (!contactElements.empty() || !mortarContactElements.empty()) {
1697 if (contactOrder > approxOrder) {
1698 Range cont_ents, cont_adj_ent;
1699 cont_ents.merge(contactSlaveFaces);
1700 cont_ents.merge(contactMasterFaces);
1701 cont_ents.merge(mortarContactMasterFaces);
1702 cont_ents.merge(mortarContactSlaveFaces);
1703
1704 CHKERR moab.get_adjacencies(cont_ents, 1, false, cont_adj_ent,
1705 moab::Interface::UNION);
1706 cont_ents.merge(cont_adj_ent);
1707 cont_ents = intersect(cont_ents, non_ref_ents);
1708 CHKERR mField.set_field_order(cont_ents, "SPATIAL_POSITION", contactOrder,
1709 verb);
1710 if (debug) {
1711 CHKERR write_file("contact_order_" + std::to_string(contactOrder) +
1712 ".vtk",
1713 cont_ents);
1714 }
1715 }
1716
1717 // add Lagrange multipliers field for contact constraints on slave faces
1718 CHKERR mField.add_field("LAMBDA_CONTACT", H1, AINSWORTH_LEGENDRE_BASE, 1,
1719 MB_TAG_SPARSE, MF_ZERO);
1720
1721 if (!contactSlaveFaces.empty())
1723 "LAMBDA_CONTACT");
1724
1725 if (!mortarContactSlaveFaces.empty())
1727 "LAMBDA_CONTACT");
1728
1729 CHKERR mField.set_field_order(0, MBTRI, "LAMBDA_CONTACT",
1731 CHKERR mField.set_field_order(0, MBEDGE, "LAMBDA_CONTACT",
1733 CHKERR mField.set_field_order(0, MBVERTEX, "LAMBDA_CONTACT", 1);
1734 }
1735
1737
1739}
1740
1742 const bool build_fields,
1743 const int verb) {
1744 const RefEntity_multiIndex *refined_ents_ptr;
1745 const FieldEntity_multiIndex *field_ents_ptr;
1747 CHKERR mField.get_ref_ents(&refined_ents_ptr);
1748 CHKERR mField.get_field_ents(&field_ents_ptr);
1749 if (mField.check_field("LAMBDA_ARC_LENGTH")) {
1750 auto vit = refined_ents_ptr->get<Ent_mi_tag>().find(arcLengthVertex);
1751 if (vit == refined_ents_ptr->get<Ent_mi_tag>().end()) {
1752 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1753 "Vertex of arc-length field not found");
1754 }
1755 *(const_cast<RefEntity *>(vit->get())->getBitRefLevelPtr()) |= bit;
1756 auto fit = field_ents_ptr->get<Unique_mi_tag>().find(
1757 FieldEntity::getLocalUniqueIdCalculate(
1758 mField.get_field_bit_number("LAMBDA_ARC_LENGTH"), arcLengthVertex));
1759 if (fit == field_ents_ptr->get<Unique_mi_tag>().end())
1760 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1761 "Arc length field entity not found");
1762
1763 MOFEM_LOG("CPWorld", Sev::noisy) << **fit;
1764 MOFEM_LOG_C("CPWorld", Sev::inform, "Arc-Length field lambda = %3.4e",
1765 (*fit)->getEntFieldData()[0]);
1767 }
1768
1769 CHKERR mField.add_field("LAMBDA_ARC_LENGTH", NOFIELD, NOBASE, 1,
1770 MB_TAG_SPARSE, MF_ZERO, verb);
1771 EntityHandle lambda_meshset = mField.get_field_meshset("LAMBDA_ARC_LENGTH");
1772 CHKERR mField.get_moab().add_entities(lambda_meshset, &arcLengthVertex, 1);
1773 // Add vertex to MoFEM database
1774 std::pair<RefEntity_multiIndex::iterator, bool> p_ent =
1775 const_cast<RefEntity_multiIndex *>(refined_ents_ptr)
1776 ->insert(boost::make_shared<RefEntity>(
1778 *(const_cast<RefEntity *>(p_ent.first->get())->getBitRefLevelPtr()) |=
1779 (bit | BitRefLevel().set(BITREFLEVEL_SIZE - 2));
1780 if (build_fields) {
1782 }
1784}
1785
1787 const bool proc_only,
1788 const bool build_fields,
1789 const int verb,
1790 const bool debug) {
1792 MeshsetsManager *meshset_manager_ptr;
1794 CHKERR mField.getInterface(meshset_manager_ptr);
1795
1796 Range level_tris;
1797 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1798 bit, BitRefLevel().set(), MBTRI, level_tris);
1799
1800 Range body_skin_proc_tris;
1801 body_skin_proc_tris = intersect(bodySkin, level_tris);
1802
1803 if (proc_only) {
1804 Range proc_tris;
1805 CHKERR moab.get_adjacencies(bitProcEnts, 2, false, proc_tris,
1806 moab::Interface::UNION);
1807 body_skin_proc_tris = intersect(body_skin_proc_tris, proc_tris);
1808 }
1809
1810 if (debug) {
1811 EntityHandle meshset;
1812 CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
1813 CHKERR mField.get_moab().add_entities(meshset, body_skin_proc_tris);
1814 std::string field_name;
1815 field_name = "out_proc_body_skin_" +
1816 boost::lexical_cast<std::string>(mField.get_comm_rank()) +
1817 ".vtk";
1818 CHKERR mField.get_moab().write_file(field_name.c_str(), "VTK", "", &meshset,
1819 1);
1820 CHKERR mField.get_moab().delete_entities(&meshset, 1);
1821 }
1822
1823 std::string field_name =
1824 "LAMBDA_SURFACE" + boost::lexical_cast<std::string>(
1825 getInterface<CPMeshCut>()->getSkinOfTheBodyId());
1826 Range dofs_nodes;
1827 CHKERR moab.get_connectivity(body_skin_proc_tris, dofs_nodes, true);
1828 if (proc_only)
1829 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(dofs_nodes,
1830 QUIET);
1831
1833 MB_TAG_SPARSE, MF_ZERO, verb);
1834
1835 Range ents_to_remove;
1837 ents_to_remove = subtract(ents_to_remove, dofs_nodes);
1838 CHKERR mField.remove_ents_from_field(field_name, ents_to_remove, verb);
1839
1840 CHKERR mField.add_ents_to_field_by_type(dofs_nodes, MBVERTEX, field_name,
1841 verb);
1842 CHKERR mField.set_field_order(dofs_nodes, field_name, 1, verb);
1843
1844 if (build_fields) {
1846 }
1847
1849}
1850
1852 const bool proc_only,
1853 const bool build_fields,
1854 const int verb,
1855 const bool debug) {
1857 auto *meshset_manager_ptr = mField.getInterface<MeshsetsManager>();
1858 auto *bit_ref_manager_ptr = mField.getInterface<BitRefManager>();
1860
1861 Range level_edges;
1862 CHKERR bit_ref_manager_ptr->getEntitiesByTypeAndRefLevel(
1863 bit, BitRefLevel().set(), MBEDGE, level_edges);
1864 Range edges_ents;
1865 if (meshset_manager_ptr->checkMeshset(
1866 getInterface<CPMeshCut>()->getEdgesBlockSet(), BLOCKSET)) {
1867 Range meshset_ents;
1868 CHKERR meshset_manager_ptr->getEntitiesByDimension(
1869 getInterface<CPMeshCut>()->getEdgesBlockSet(), BLOCKSET, 1,
1870 meshset_ents, true);
1871 edges_ents = intersect(meshset_ents, level_edges);
1872 }
1873
1874 if (proc_only) {
1875 Range proc_edges;
1876 CHKERR moab.get_adjacencies(bitProcEnts, 1, false, proc_edges,
1877 moab::Interface::UNION);
1878 edges_ents = intersect(edges_ents, proc_edges);
1879 }
1880
1881 // Get nodes on other side faces
1882 Range other_side_crack_faces_nodes;
1883 rval = mField.get_moab().get_connectivity(otherSideCrackFaces,
1884 other_side_crack_faces_nodes, true);
1885 other_side_crack_faces_nodes =
1886 subtract(other_side_crack_faces_nodes, crackFrontNodes);
1887
1888 Range dofs_nodes;
1889 CHKERR moab.get_connectivity(edges_ents, dofs_nodes, true);
1890 dofs_nodes = subtract(dofs_nodes, other_side_crack_faces_nodes);
1891 if (proc_only)
1892 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(dofs_nodes,
1893 QUIET);
1894
1895 std::string field_name = "LAMBDA_EDGE";
1897 MB_TAG_SPARSE, MF_ZERO, verb);
1898
1899 Range ents_to_remove;
1901 ents_to_remove = subtract(ents_to_remove, dofs_nodes);
1902 CHKERR mField.remove_ents_from_field(field_name, ents_to_remove, verb);
1903
1904 CHKERR mField.add_ents_to_field_by_type(dofs_nodes, MBVERTEX, field_name,
1905 verb);
1906 CHKERR mField.set_field_order(dofs_nodes, field_name, 1, verb);
1907
1908 auto get_fields_multi_index = [this]() {
1909 const Field_multiIndex *fields_ptr;
1910 CHKERR mField.get_fields(&fields_ptr);
1911 return *fields_ptr;
1912 };
1913 for (auto field : get_fields_multi_index()) {
1914 if (field->getName().compare(0, 14, "LAMBDA_SURFACE") == 0) {
1915 CHKERR mField.remove_ents_from_field(field->getName(), dofs_nodes, verb);
1916 }
1917 }
1918
1919 if (build_fields) {
1921 }
1922
1924}
1925
1927 const BitRefLevel bit, const bool proc_only, const bool build_fields,
1928 const int verb, const bool debug) {
1930 MeshsetsManager *meshset_manager_ptr;
1932
1933 CHKERR mField.getInterface(meshset_manager_ptr);
1934 Range crack_surface_ents = oneSideCrackFaces;
1935 Range level_tris;
1936 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1937 bit, BitRefLevel().set(), MBTRI, level_tris);
1938 crack_surface_ents = intersect(crack_surface_ents, level_tris);
1939
1940 if (proc_only) {
1941 Range proc_tris;
1942 CHKERR moab.get_adjacencies(bitProcEnts, 2, false, proc_tris,
1943 moab::Interface::UNION);
1944 crack_surface_ents = intersect(crack_surface_ents, proc_tris);
1945 }
1946
1947 if (debug) {
1948 EntityHandle meshset;
1949 CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
1950 CHKERR mField.get_moab().add_entities(meshset, crack_surface_ents);
1951 std::string field_name;
1952 field_name = "out_proc_crack_" +
1953 boost::lexical_cast<std::string>(mField.get_comm_rank()) +
1954 ".vtk";
1955 CHKERR mField.get_moab().write_file(field_name.c_str(), "VTK", "", &meshset,
1956 1);
1957 CHKERR mField.get_moab().delete_entities(&meshset, 1);
1958 }
1959
1960 std::string field_name =
1961 "LAMBDA_SURFACE" + boost::lexical_cast<std::string>(
1962 getInterface<CPMeshCut>()->getCrackSurfaceId());
1963 Range dofs_nodes;
1964 CHKERR moab.get_connectivity(crack_surface_ents, dofs_nodes, true);
1965 dofs_nodes = subtract(dofs_nodes, crackFrontNodes);
1966 if (proc_only)
1967 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(dofs_nodes,
1968 QUIET);
1969
1971 MB_TAG_SPARSE, MF_ZERO, verb);
1972 Range ents_to_remove;
1974 ents_to_remove = subtract(ents_to_remove, dofs_nodes);
1975 CHKERR mField.remove_ents_from_field(field_name, ents_to_remove, verb);
1976
1977 CHKERR mField.add_ents_to_field_by_type(dofs_nodes, MBVERTEX, field_name,
1978 verb);
1979 CHKERR mField.set_field_order(dofs_nodes, field_name, 1, verb);
1980
1981 if (build_fields) {
1983 }
1984
1986}
1987
1989 const BitRefLevel bit_spatial, const BitRefLevel bit_material,
1990 const bool proc_only, const bool build_fields, const int verb,
1991 const bool debug) {
1993
1994 auto add_both_sides_field = [&](const std::string lambda_field_name,
1995 Range master_nodes, Range nodes_to_remove,
1996 bool add_ho) {
1998 // This field set constrains that both side of crack surface have the same
1999 // node positions
2000 CHKERR mField.add_field(lambda_field_name, H1, AINSWORTH_LEGENDRE_BASE, 3,
2001 MB_TAG_SPARSE, MF_ZERO);
2002
2003 if (debug) {
2004 std::ostringstream file;
2005 file << lambda_field_name << ".vtk";
2006 EntityHandle meshset;
2007 CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
2008 CHKERR mField.get_moab().add_entities(meshset, master_nodes);
2009 CHKERR mField.get_moab().write_file(file.str().c_str(), "VTK", "",
2010 &meshset, 1);
2011 CHKERR mField.get_moab().delete_entities(&meshset, 1);
2012 }
2013
2014 // remove old ents from field which are no longer part of the field
2015 Range ents_to_remove;
2016 CHKERR mField.get_field_entities_by_handle(lambda_field_name,
2017 ents_to_remove);
2018 ents_to_remove = subtract(ents_to_remove, master_nodes);
2019 CHKERR mField.remove_ents_from_field(lambda_field_name, ents_to_remove,
2020 verb);
2021
2022 CHKERR mField.add_ents_to_field_by_type(master_nodes, MBVERTEX,
2023 lambda_field_name);
2024 if (approxOrder > 1 && add_ho) {
2025 CHKERR mField.add_ents_to_field_by_type(master_nodes, MBEDGE,
2026 lambda_field_name);
2027 CHKERR mField.add_ents_to_field_by_type(master_nodes, MBTRI,
2028 lambda_field_name);
2029 }
2030 CHKERR mField.set_field_order(master_nodes.subset_by_dimension(0),
2031 lambda_field_name, 1);
2032 if (approxOrder > 1 && add_ho) {
2033 CHKERR mField.set_field_order(master_nodes.subset_by_dimension(1),
2034 lambda_field_name, approxOrder);
2035 CHKERR mField.set_field_order(master_nodes.subset_by_dimension(2),
2036 lambda_field_name, approxOrder);
2037 }
2038
2039 // Since we state constrains that mesh nodes, in material space, have to
2040 // move the same on the top and bottom side of crack surface, we have to
2041 // remove constraint from the all surface constrains fields, except
2042 // "LAMBDA_SURFACE", otherwise system will be over-constrained.
2043 const Field_multiIndex *fields_ptr;
2044 CHKERR mField.get_fields(&fields_ptr);
2045 for (auto const &field : (*fields_ptr)) {
2046 const std::string field_name = field->getName();
2047 if (field_name.compare(0, 14, "LAMBDA_SURFACE") == 0) {
2048 CHKERR mField.remove_ents_from_field(field_name, nodes_to_remove, verb);
2049 }
2050 if (field_name.compare(0, 11, "LAMBDA_EDGE") == 0) {
2051 CHKERR mField.remove_ents_from_field(field_name, nodes_to_remove, verb);
2052 }
2053 }
2055 };
2056
2057 auto write_file = [&](const std::string file_name, Range ents) {
2059 EntityHandle meshset;
2060 CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
2061 CHKERR mField.get_moab().add_entities(meshset, ents);
2062 if (mField.get_comm_rank() == 0) {
2063 CHKERR mField.get_moab().write_file(file_name.c_str(), "VTK", "",
2064 &meshset, 1);
2065 }
2066 CHKERR mField.get_moab().delete_entities(&meshset, 1);
2068 };
2069
2070 auto get_prisms_level_nodes = [&](auto bit) {
2071 Range prisms_level;
2072 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2073 bit, BitRefLevel().set(), MBPRISM, prisms_level);
2074 if (proc_only) {
2075 prisms_level = intersect(bitProcEnts, prisms_level);
2076 }
2077 Range prisms_level_nodes;
2078 CHKERR mField.get_moab().get_connectivity(prisms_level, prisms_level_nodes,
2079 true);
2080 return prisms_level_nodes;
2081 };
2082
2083 auto get_other_side_crack_faces_nodes = [&](auto prisms_level_nodes) {
2084 Range other_side_crack_faces_nodes;
2085 CHKERR mField.get_moab().get_connectivity(
2086 otherSideCrackFaces, other_side_crack_faces_nodes, true);
2087 other_side_crack_faces_nodes =
2088 intersect(other_side_crack_faces_nodes, prisms_level_nodes);
2089 other_side_crack_faces_nodes =
2090 subtract(other_side_crack_faces_nodes, crackFrontNodes);
2091
2092 return other_side_crack_faces_nodes;
2093 };
2094
2095 auto prisms_level_nodes = get_prisms_level_nodes(bit_material);
2096 auto other_side_crack_faces_nodes =
2097 get_other_side_crack_faces_nodes(prisms_level_nodes);
2098
2099 CHKERR add_both_sides_field("LAMBDA_BOTH_SIDES", other_side_crack_faces_nodes,
2100 other_side_crack_faces_nodes, false);
2101
2103 Range other_side_crack_faces_edges;
2104 if (approxOrder > 1) {
2105 CHKERR mField.get_moab().get_adjacencies(otherSideCrackFaces, 1, false,
2106 other_side_crack_faces_edges,
2107 moab::Interface::UNION);
2108 other_side_crack_faces_edges =
2109 subtract(other_side_crack_faces_edges, crackFront);
2110 other_side_crack_faces_edges.merge(otherSideCrackFaces);
2111 }
2112 other_side_crack_faces_edges.merge(
2113 get_other_side_crack_faces_nodes(get_prisms_level_nodes(bit_spatial)));
2114 CHKERR add_both_sides_field("LAMBDA_CLOSE_CRACK",
2115 other_side_crack_faces_edges, Range(), true);
2116 }
2117
2118 if (!contactElements.empty() && !ignoreContact && !fixContactNodes) {
2120 CHKERR mField.get_moab().get_connectivity(
2123 intersect(contactBothSidesMasterNodes, prisms_level_nodes);
2125 subtract(contactBothSidesMasterNodes, other_side_crack_faces_nodes);
2126 Range contact_both_sides_slave_nodes;
2127 for (Range::iterator nit = contactBothSidesMasterNodes.begin();
2128 nit != contactBothSidesMasterNodes.end(); nit++) {
2129 Range contact_prisms;
2130 CHKERR mField.get_moab().get_adjacencies(
2131 &*nit, 1, 3, false, contact_prisms, moab::Interface::UNION);
2132 contact_prisms = intersect(contact_prisms, contactElements);
2133 EntityHandle prism = contact_prisms.front();
2134 const EntityHandle *conn;
2135 int side_number, other_side_number, sense, offset, number_nodes = 0;
2136 CHKERR mField.get_moab().get_connectivity(prism, conn, number_nodes);
2137 CHKERR mField.get_moab().side_number(prism, *nit, side_number, sense,
2138 offset);
2139 if (side_number < 3)
2140 contact_both_sides_slave_nodes.insert(conn[side_number + 3]);
2141 else
2142 contact_both_sides_slave_nodes.insert(conn[side_number - 3]);
2143 }
2144 contact_both_sides_slave_nodes =
2145 intersect(contact_both_sides_slave_nodes, prisms_level_nodes);
2146
2147 if (debug) {
2148 CHKERR write_file("contact_both_sides_slave_nodes.vtk",
2149 contact_both_sides_slave_nodes);
2150 }
2151
2152 CHKERR add_both_sides_field("LAMBDA_BOTH_SIDES_CONTACT",
2154 contact_both_sides_slave_nodes, false);
2155 }
2156
2158}
2159
2161 const bool build_fields,
2162 const int verb,
2163 const bool debug) {
2166
2167 Range level_edges;
2168 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2169 bit, BitRefLevel().set(), MBEDGE, level_edges);
2170 Range crack_front_edges = intersect(crackFront, level_edges);
2171 Range crack_front_nodes;
2172 CHKERR moab.get_connectivity(crack_front_edges, crack_front_nodes, true);
2173
2174 CHKERR mField.add_field("LAMBDA_CRACKFRONT_AREA", H1, AINSWORTH_LEGENDRE_BASE,
2175 1, MB_TAG_SPARSE, MF_ZERO, verb);
2176 // Add dofs to LAMBDA_CRACKFRONT_AREA field
2177 {
2178 Range ents_to_remove;
2179 CHKERR mField.get_field_entities_by_handle("LAMBDA_CRACKFRONT_AREA",
2180 ents_to_remove);
2181 ents_to_remove = subtract(ents_to_remove, crack_front_nodes);
2182 CHKERR mField.remove_ents_from_field("LAMBDA_CRACKFRONT_AREA",
2183 ents_to_remove, verb);
2184 CHKERR mField.add_ents_to_field_by_type(crack_front_nodes, MBVERTEX,
2185 "LAMBDA_CRACKFRONT_AREA", verb);
2186 CHKERR mField.set_field_order(crack_front_nodes, "LAMBDA_CRACKFRONT_AREA",
2187 1, verb);
2188 }
2189
2190 CHKERR mField.add_field("LAMBDA_CRACKFRONT_AREA_TANGENT", H1,
2191 AINSWORTH_LEGENDRE_BASE, 1, MB_TAG_SPARSE, MF_ZERO,
2192 verb);
2193 // Add dofs to LAMBDA_CRACKFRONT_AREA_TANGENT field
2194 {
2195 // Remove dofs on crack front ends
2196 Skinner skin(&mField.get_moab());
2197 Range skin_crack_front;
2198 rval = skin.find_skin(0, crack_front_edges, false, skin_crack_front);
2199 Range field_ents = subtract(crack_front_nodes, skin_crack_front);
2200 // remove old ents from field which are no longer part of the field
2201 Range ents_to_remove;
2202 CHKERR mField.get_field_entities_by_handle("LAMBDA_CRACKFRONT_AREA_TANGENT",
2203 ents_to_remove);
2204 ents_to_remove = subtract(ents_to_remove, field_ents);
2205 CHKERR mField.remove_ents_from_field("LAMBDA_CRACKFRONT_AREA_TANGENT",
2206 ents_to_remove, verb);
2208 field_ents, MBVERTEX, "LAMBDA_CRACKFRONT_AREA_TANGENT", verb);
2209 CHKERR mField.set_field_order(field_ents, "LAMBDA_CRACKFRONT_AREA_TANGENT",
2210 1, verb);
2211 }
2212
2213 if (build_fields) {
2215 }
2216
2218}
2219
2221 const BitRefLevel bit1, const BitRefLevel mask1, const BitRefLevel bit2,
2222 const BitRefLevel mask2, const bool add_forces, const bool proc_only,
2223 const int verb) {
2225
2226 Range tets_level1;
2227 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2228 bit1, mask1, MBTET, tets_level1);
2229 if (proc_only) {
2230 tets_level1 = intersect(bitProcEnts, tets_level1);
2231 }
2232 Range tets_not_level2;
2233 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2234 bit2, mask2, MBTET, tets_not_level2);
2235 tets_not_level2 = subtract(tets_level1, tets_not_level2);
2236
2237 // Add finite elements
2238 CHKERR mField.add_finite_element("ELASTIC", MF_ZERO, verb);
2240 "SPATIAL_POSITION");
2242 "SPATIAL_POSITION");
2244 "SPATIAL_POSITION");
2246 "MESH_NODE_POSITIONS");
2247 // Add ents and build finite elements. Only build elements which has been
2248 // added.
2249 {
2250 Range current_ents_with_fe;
2252 current_ents_with_fe);
2253 Range ents_to_remove;
2254 ents_to_remove = subtract(current_ents_with_fe, tets_level1);
2255
2256 CHKERR mField.remove_ents_from_finite_element("ELASTIC", ents_to_remove);
2258 "ELASTIC");
2259 CHKERR mField.build_finite_elements("ELASTIC", &tets_level1, VERBOSE);
2260 }
2261
2262 // Add finite elements not on crack front
2263 CHKERR mField.add_finite_element("ELASTIC_NOT_ALE", MF_ZERO, verb);
2265 "SPATIAL_POSITION");
2267 "SPATIAL_POSITION");
2269 "SPATIAL_POSITION");
2271 "MESH_NODE_POSITIONS");
2272
2273 // Add ents and build finite elements
2274 {
2275 Range current_ents_with_fe;
2277 current_ents_with_fe);
2278 Range ents_to_remove;
2279 ents_to_remove = subtract(current_ents_with_fe, tets_not_level2);
2280
2282 ents_to_remove);
2283 CHKERR mField.add_ents_to_finite_element_by_type(tets_not_level2, MBTET,
2284 "ELASTIC_NOT_ALE");
2285 CHKERR mField.build_finite_elements("ELASTIC_NOT_ALE", &tets_not_level2,
2286 verb);
2287 }
2288
2289 // Add skin elements for faster postprocessing
2290 {
2291 CHKERR mField.add_finite_element("SKIN", MF_ZERO, verb);
2292 Skinner skin(&mField.get_moab());
2293 Range proc_tris, skin_tris, tets_level;
2294 CHKERR mField.get_moab().get_adjacencies(bitProcEnts, 2, false, proc_tris,
2295 moab::Interface::UNION);
2296
2297 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2298 bit1, mask1, MBTET, tets_level);
2299
2300 CHKERR skin.find_skin(0, tets_level, false, skin_tris);
2301 skin_tris = intersect(skin_tris, proc_tris);
2302
2303 CHKERR mField.add_ents_to_finite_element_by_type(skin_tris, MBTRI, "SKIN");
2304
2306 "SPATIAL_POSITION");
2308 "SPATIAL_POSITION");
2310 "SPATIAL_POSITION");
2312 "MESH_NODE_POSITIONS");
2315 "SKIN", "EIGEN_SPATIAL_POSITIONS");
2317 }
2318 // Add spring boundary condition applied on surfaces.
2319 // This is only declaration not implementation.
2320 CHKERR MetaSpringBC::addSpringElements(mField, "SPATIAL_POSITION",
2321 "MESH_NODE_POSITIONS");
2323
2324 auto cn_value_ptr = boost::make_shared<double>(cnValue);
2325
2326 contactProblem = boost::shared_ptr<SimpleContactProblem>(
2327 new SimpleContactProblem(mField, cn_value_ptr));
2328 // add fields to the global matrix by adding the element
2329
2330 CHKERR contactProblem->addContactElement(
2331 "CONTACT", "SPATIAL_POSITION", "LAMBDA_CONTACT", contactElements,
2332 solveEigenStressProblem, "EIGEN_SPATIAL_POSITIONS");
2333
2335
2337 boost::shared_ptr<MortarContactProblem>(new MortarContactProblem(
2338 mField, contactSearchMultiIndexPtr, cn_value_ptr));
2339
2340 CHKERR mortarContactProblemPtr->addMortarContactElement(
2341 "MORTAR_CONTACT", "SPATIAL_POSITION", "LAMBDA_CONTACT",
2342 mortarContactElements, "MESH_NODE_POSITIONS", solveEigenStressProblem,
2343 "EIGEN_SPATIAL_POSITIONS");
2344
2345 CHKERR mField.build_finite_elements("MORTAR_CONTACT");
2346
2347 Range all_slave_tris;
2348 all_slave_tris.merge(mortarContactSlaveFaces);
2349 all_slave_tris.merge(contactSlaveFaces);
2350
2351 CHKERR contactProblem->addPostProcContactElement(
2352 "CONTACT_POST_PROC", "SPATIAL_POSITION", "LAMBDA_CONTACT",
2353 "MESH_NODE_POSITIONS", all_slave_tris);
2354 CHKERR mField.build_finite_elements("CONTACT_POST_PROC");
2355
2357}
2358
2360 const int verb) {
2361 const EntFiniteElement_multiIndex *fe_ent_ptr;
2363
2365 auto &fe_ent_by_name = fe_ent_ptr->get<FiniteElement_name_mi_tag>();
2366
2367 if (mField.check_finite_element("ARC_LENGTH")) {
2368 Range arc_fe_meshsets;
2369 for (auto fit = fe_ent_by_name.lower_bound("ARC_LENGTH");
2370 fit != fe_ent_by_name.upper_bound("ARC_LENGTH"); ++fit) {
2371 arc_fe_meshsets.insert(fit->get()->getEnt());
2372 if (auto ptr = fit->get()->getPartProcPtr())
2373 *ptr = 0;
2374 }
2375 CHKERR mField.getInterface<BitRefManager>()->addBitRefLevel(arc_fe_meshsets,
2376 bit);
2377 CHKERR mField.build_finite_elements("ARC_LENGTH", NULL, verb);
2379 }
2380
2381 CHKERR mField.add_finite_element("ARC_LENGTH", MF_ZERO, verb);
2383 "LAMBDA_ARC_LENGTH");
2385 "LAMBDA_ARC_LENGTH");
2386 // elem data
2388 "LAMBDA_ARC_LENGTH");
2389
2390 // finally add created meshset to the ARC_LENGTH finite element
2391 EntityHandle meshset_fe_arc_length;
2392 {
2393 CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset_fe_arc_length);
2394 CHKERR mField.get_moab().add_entities(meshset_fe_arc_length,
2395 &arcLengthVertex, 1);
2396 CHKERR mField.getInterface<BitRefManager>()->setBitLevelToMeshset(
2397 meshset_fe_arc_length, bit | BitRefLevel().set(BITREFLEVEL_SIZE - 2));
2398 }
2400 "ARC_LENGTH", false);
2401 CHKERR mField.build_finite_elements("ARC_LENGTH", NULL, verb);
2402 for (auto fit = fe_ent_by_name.lower_bound("ARC_LENGTH");
2403 fit != fe_ent_by_name.upper_bound("ARC_LENGTH"); ++fit) {
2404 if (auto ptr = fit->get()->getPartProcPtr())
2405 *ptr = 0;
2406 }
2408}
2409
2411 const BitRefLevel mask,
2412 const bool proc_only) {
2415
2416 Range tets_level;
2417 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2418 bit, mask, MBTET, tets_level);
2419 if (proc_only) {
2420 tets_level = intersect(bitProcEnts, tets_level);
2421 }
2422
2423 // Set approx order
2424 Range dofs_ents;
2425 dofs_ents.merge(tets_level);
2426 for (int dd = 1; dd != 3; dd++) {
2427 CHKERR moab.get_adjacencies(tets_level, dd, false, dofs_ents,
2428 moab::Interface::UNION);
2429 }
2430 Range dofs_nodes;
2431 rval = moab.get_connectivity(tets_level, dofs_nodes, true);
2433 Range lower_dim_ents = unite(dofs_ents, dofs_nodes);
2434
2435#ifdef __ANALITICAL_TRACTION__
2436 {
2437 CHKERR mField.add_finite_element("ANALITICAL_TRACTION", MF_ZERO);
2439 "SPATIAL_POSITION");
2441 "SPATIAL_POSITION");
2443 "SPATIAL_POSITION");
2445 "MESH_NODE_POSITIONS");
2446 MeshsetsManager *meshset_manager_ptr;
2447 CHKERR mField.getInterface(meshset_manager_ptr);
2448 if (meshset_manager_ptr->checkMeshset("ANALITICAL_TRACTION")) {
2449 const CubitMeshSets *cubit_meshset_ptr;
2450 meshset_manager_ptr->getCubitMeshsetPtr("ANALITICAL_TRACTION",
2451 &cubit_meshset_ptr);
2452 Range tris;
2453 CHKERR meshset_manager_ptr->getEntitiesByDimension(
2454 cubit_meshset_ptr->getMeshsetId(), BLOCKSET, 2, tris, true);
2455 tris = intersect(tris, lower_dim_ents);
2457 "ANALITICAL_TRACTION");
2458 }
2459 CHKERR mField.build_finite_elements("ANALITICAL_TRACTION", &lower_dim_ents);
2460 }
2461#endif //__ANALITICAL_TRACTION__
2462
2465 "SPATIAL_POSITION");
2467 "SPATIAL_POSITION");
2469 "SPATIAL_POSITION");
2471 "MESH_NODE_POSITIONS");
2473 it)) {
2474 Range tris;
2475 CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
2476 true);
2477 Range edges;
2478 CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBEDGE, edges,
2479 true);
2480 Range tris_edges;
2481 CHKERR mField.get_moab().get_adjacencies(tris, 1, false, tris_edges,
2482 moab::Interface::UNION);
2483 Range tris_nodes;
2484 CHKERR mField.get_moab().get_connectivity(tris, tris_nodes);
2485 Range edges_nodes;
2486 CHKERR mField.get_moab().get_connectivity(edges, edges_nodes);
2487 Range nodes;
2488 CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBVERTEX, nodes,
2489 true);
2490 nodes = subtract(nodes, tris_nodes);
2491 nodes = subtract(nodes, edges_nodes);
2492 nodes = intersect(lower_dim_ents, nodes);
2493 edges = subtract(edges, tris_edges);
2494 edges = intersect(lower_dim_ents, edges);
2495 tris = intersect(lower_dim_ents, tris);
2496 CHKERR mField.add_ents_to_finite_element_by_type(tris, MBTRI, "FORCE_FE");
2497 CHKERR mField.add_ents_to_finite_element_by_type(edges, MBEDGE, "FORCE_FE");
2499 "FORCE_FE");
2500 }
2501
2502 // Reading forces from BLOCKSET (for bone meshes)
2503 const string block_set_force_name("FORCE");
2504 // search for block named FORCE and add its attributes to FORCE_FE element
2506 if (it->getName().compare(0, block_set_force_name.length(),
2507 block_set_force_name) == 0) {
2508 std::vector<double> mydata;
2509 CHKERR it->getAttributes(mydata);
2510
2511 VectorDouble force(mydata.size());
2512 for (unsigned int ii = 0; ii < mydata.size(); ii++) {
2513 force[ii] = mydata[ii];
2514 }
2515 if (force.size() != 3) {
2516 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2517 "Force vector not given");
2518 }
2519 Range tris;
2520 CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
2521 true);
2522 CHKERR mField.add_ents_to_finite_element_by_type(tris, MBTRI, "FORCE_FE");
2523 }
2524 }
2525
2526 // CHKERR MetaNodalForces::addElement(mField,"SPATIAL_POSITION");
2527 // CHKERR MetaEdgeForces::addElement(mField,"SPATIAL_POSITION");
2528 CHKERR mField.build_finite_elements("FORCE_FE", &lower_dim_ents);
2529 CHKERR mField.add_finite_element("PRESSURE_FE", MF_ZERO);
2531 "SPATIAL_POSITION");
2533 "SPATIAL_POSITION");
2535 "SPATIAL_POSITION");
2537 "MESH_NODE_POSITIONS");
2539 it)) {
2540 Range tris;
2541 CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
2542 true);
2543 tris = intersect(lower_dim_ents, tris);
2545 "PRESSURE_FE");
2546 }
2547 const string block_set_linear_pressure_name("LINEAR_PRESSURE");
2549 if (it->getName().compare(0, block_set_linear_pressure_name.length(),
2550 block_set_linear_pressure_name) == 0) {
2551 Range tris;
2552 CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
2553 true);
2554 tris = intersect(lower_dim_ents, tris);
2556 "PRESSURE_FE");
2557 }
2558 }
2559
2560 // search for block named PRESSURE and add its attributes to PRESSURE_FE
2561 // element (for bone meshes)
2562 const string block_set_pressure_name("PRESSURE");
2564 if (it->getName().compare(0, block_set_pressure_name.length(),
2565 block_set_pressure_name) == 0) {
2566 std::vector<double> mydata;
2567 CHKERR it->getAttributes(mydata);
2568 VectorDouble pressure(mydata.size());
2569 for (unsigned int ii = 0; ii < mydata.size(); ii++) {
2570 pressure[ii] = mydata[ii];
2571 }
2572 if (pressure.empty()) {
2573 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2574 "Pressure not given");
2575 }
2576 Range tris;
2577 CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
2578 true);
2579
2580 CHKERR
2581 mField.add_ents_to_finite_element_by_type(tris, MBTRI, "PRESSURE_FE");
2582 }
2583 }
2584
2585 CHKERR mField.build_finite_elements("PRESSURE_FE", &lower_dim_ents);
2586
2587#ifdef __ANALITICAL_DISPLACEMENT__
2588 {
2589 MeshsetsManager *meshset_manager_ptr;
2590 CHKERR mField.getInterface(meshset_manager_ptr);
2591 if (meshset_manager_ptr->checkMeshset("ANALITICAL_DISP")) {
2592 const CubitMeshSets *cubit_meshset_ptr;
2593 meshset_manager_ptr->getCubitMeshsetPtr("ANALITICAL_DISP",
2594 &cubit_meshset_ptr);
2595 Range tris;
2596 CHKERR meshset_manager_ptr->getEntitiesByDimension(
2597 cubit_meshset_ptr->getMeshsetId(), BLOCKSET, 2, tris, true);
2598 if (!tris.empty()) {
2601 "SPATIAL_POSITION");
2603 "SPATIAL_POSITION");
2605 "SPATIAL_POSITION");
2607 "BC_FE", "MESH_NODE_POSITIONS");
2608 tris = intersect(lower_dim_ents, tris);
2609 CHKERR mField.add_ents_to_finite_element_by_type(tris, MBTRI, "BC_FE");
2610 CHKERR mField.build_finite_elements("BC_FE", &tris);
2611 }
2612 }
2613 }
2614#endif //__ANALITICAL_DISPLACEMENT__
2615
2617}
2618
2620 const BitRefLevel mask,
2621 const bool proc_only) {
2622
2625
2626 Range tets_level;
2627 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2628 bit, mask, MBTET, tets_level);
2629 if (proc_only) {
2630 tets_level = intersect(bitProcEnts, tets_level);
2631 }
2632
2633 Range dofs_ents;
2634 dofs_ents.merge(tets_level);
2635 for (int dd = 1; dd != 3; dd++) {
2636 CHKERR moab.get_adjacencies(tets_level, dd, false, dofs_ents,
2637 moab::Interface::UNION);
2638 }
2639 Range dofs_nodes;
2640 rval = moab.get_connectivity(tets_level, dofs_nodes, true);
2642 Range lower_dim_ents = unite(dofs_ents, dofs_nodes);
2643
2644 CHKERR mField.add_finite_element("PRESSURE_ALE", MF_ZERO);
2646 "SPATIAL_POSITION");
2648 "SPATIAL_POSITION");
2650 "SPATIAL_POSITION");
2652 "MESH_NODE_POSITIONS");
2654 "MESH_NODE_POSITIONS");
2656 "MESH_NODE_POSITIONS");
2657
2658 Range ents_to_add;
2660 it)) {
2661 Range tris;
2662 CHKERR mField.get_moab().get_entities_by_type(it->meshset, MBTRI, tris,
2663 true);
2664 tris = intersect(lower_dim_ents, tris);
2665 ents_to_add.merge(tris);
2666 }
2667
2668 Range current_ents_with_fe;
2670 current_ents_with_fe);
2671 Range ents_to_remove;
2672 ents_to_remove = subtract(current_ents_with_fe, ents_to_add);
2673 CHKERR mField.remove_ents_from_finite_element("PRESSURE_ALE", ents_to_remove);
2675 "PRESSURE_ALE");
2676
2677 CHKERR mField.build_finite_elements("PRESSURE_ALE", &lower_dim_ents);
2679}
2680
2682 const BitRefLevel mask,
2683 const bool proc_only) {
2684
2687
2688 Range tets_level;
2689 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2690 bit, mask, MBTET, tets_level);
2691 if (proc_only) {
2692 tets_level = intersect(bitProcEnts, tets_level);
2693 }
2694
2695 Range dofs_ents;
2696 dofs_ents.merge(tets_level);
2697 for (int dd = 1; dd != 3; dd++) {
2698 CHKERR moab.get_adjacencies(tets_level, dd, false, dofs_ents,
2699 moab::Interface::UNION);
2700 }
2701 Range dofs_nodes;
2702 rval = moab.get_connectivity(tets_level, dofs_nodes, true);
2704 Range lower_dim_ents = unite(dofs_ents, dofs_nodes);
2705
2706 CHKERR mField.add_finite_element("SPRING_ALE", MF_ZERO);
2708 "SPATIAL_POSITION");
2710 "SPATIAL_POSITION");
2712 "SPATIAL_POSITION");
2714 "MESH_NODE_POSITIONS");
2716 "MESH_NODE_POSITIONS");
2718 "MESH_NODE_POSITIONS");
2719
2720 Range ents_to_add;
2722 if (bit->getName().compare(0, 9, "SPRING_BC") == 0) {
2723 Range tris;
2724 CHKERR mField.get_moab().get_entities_by_type(bit->meshset, MBTRI, tris,
2725 true);
2726 tris = intersect(lower_dim_ents, tris);
2727 ents_to_add.merge(tris);
2728 }
2729 }
2730
2731 Range current_ents_with_fe;
2733 current_ents_with_fe);
2734 Range ents_to_remove;
2735 ents_to_remove = subtract(current_ents_with_fe, ents_to_add);
2736 CHKERR mField.remove_ents_from_finite_element("SPRING_ALE", ents_to_remove);
2738 "SPRING_ALE");
2739
2740 CHKERR mField.build_finite_elements("SPRING_ALE", &lower_dim_ents);
2742}
2743
2745 const BitRefLevel bit, const BitRefLevel mask, const bool proc_only) {
2746
2749
2750 Range prisms_level;
2751 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2752 bit, mask, MBPRISM, prisms_level);
2753
2754 if (proc_only) {
2755 prisms_level = intersect(bitProcEnts, prisms_level);
2756 }
2757
2758 Range contact_prisms = intersect(prisms_level, contactElements);
2759
2760 CHKERR contactProblem->addContactElementALE(
2761 "SIMPLE_CONTACT_ALE", "SPATIAL_POSITION", "LAMBDA_CONTACT",
2762 "MESH_NODE_POSITIONS", contact_prisms, solveEigenStressProblem,
2763 "EIGEN_SPATIAL_POSITIONS");
2764
2765 CHKERR mField.build_finite_elements("SIMPLE_CONTACT_ALE", &contact_prisms);
2766
2767 Range face_on_prism;
2768 CHKERR moab.get_adjacencies(contact_prisms, 2, false, face_on_prism,
2769 moab::Interface::UNION);
2770
2771 Range tris_on_prism = face_on_prism.subset_by_type(MBTRI);
2772
2773 Range check_tris =
2774 intersect(tris_on_prism, unite(contactSlaveFaces, contactMasterFaces));
2775
2776 Range adj_vols_to_prisms;
2777 CHKERR moab.get_adjacencies(check_tris, 3, false, adj_vols_to_prisms,
2778 moab::Interface::UNION);
2779
2780 Range tet_level;
2781 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2782 bit, mask, MBTET, tet_level);
2783
2784 Range contact_tets_from_vols = adj_vols_to_prisms.subset_by_type(MBTET);
2785
2786 Range contact_tets = intersect(tet_level, contact_tets_from_vols);
2787
2788 CHKERR mField.add_finite_element("MAT_CONTACT", MF_ZERO);
2790 "SPATIAL_POSITION");
2792 "SPATIAL_POSITION");
2794 "MESH_NODE_POSITIONS");
2796 "MESH_NODE_POSITIONS");
2798 "SPATIAL_POSITION");
2800 "MESH_NODE_POSITIONS");
2801
2802 Range current_ale_tets;
2804 current_ale_tets);
2805 Range ale_tets_to_remove;
2806 ale_tets_to_remove = subtract(current_ale_tets, contact_tets);
2808 ale_tets_to_remove);
2809
2811 "MAT_CONTACT");
2812 CHKERR mField.build_finite_elements("MAT_CONTACT", &contact_tets);
2813
2815}
2816
2818 const BitRefLevel mask,
2819 const bool proc_only,
2820 const bool verb) {
2822
2823 Range tets_level;
2824 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2825 bit, mask, MBTET, tets_level);
2826 if (proc_only) {
2827 tets_level = intersect(bitProcEnts, tets_level);
2828 }
2829
2830 // Add finite elements
2831 CHKERR mField.add_finite_element("MATERIAL", MF_ZERO, verb);
2833 "SPATIAL_POSITION");
2835 "SPATIAL_POSITION");
2837 "MESH_NODE_POSITIONS");
2839 "MESH_NODE_POSITIONS");
2841 "SPATIAL_POSITION");
2843 "MESH_NODE_POSITIONS");
2844
2845 {
2846 Range current_ents_with_fe;
2848 current_ents_with_fe);
2849 Range ents_to_remove;
2850 ents_to_remove = subtract(current_ents_with_fe, tets_level);
2851 CHKERR mField.remove_ents_from_finite_element("MATERIAL", ents_to_remove);
2853 "MATERIAL");
2854 // Range ents_to_build = subtract(tets_level, current_ents_with_fe);
2855 // CHKERR mField.build_finite_elements("MATERIAL", &ents_to_build);
2856 CHKERR mField.build_finite_elements("MATERIAL", &tets_level, verb);
2857 }
2858
2860}
2861
2863 const BitRefLevel mask,
2864 const bool proc_only,
2865 const bool verb) {
2868
2869 CHKERR mField.add_finite_element("GRIFFITH_FORCE_ELEMENT", MF_ZERO, verb);
2870 CHKERR mField.modify_finite_element_add_field_row("GRIFFITH_FORCE_ELEMENT",
2871 "MESH_NODE_POSITIONS");
2872 CHKERR mField.modify_finite_element_add_field_col("GRIFFITH_FORCE_ELEMENT",
2873 "MESH_NODE_POSITIONS");
2874 CHKERR mField.modify_finite_element_add_field_data("GRIFFITH_FORCE_ELEMENT",
2875 "MESH_NODE_POSITIONS");
2876 Range crack_front_adj_tris;
2877 CHKERR mField.get_moab().get_adjacencies(
2878 crackFrontNodes, 2, false, crack_front_adj_tris, moab::Interface::UNION);
2879 crack_front_adj_tris = intersect(crack_front_adj_tris, crackFaces);
2880 Range tris_level;
2881 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2882 bit, mask, MBTRI, tris_level);
2883 crack_front_adj_tris = intersect(crack_front_adj_tris, tris_level);
2884 crackFrontNodesTris = crack_front_adj_tris;
2885 Range ents_to_remove;
2886 CHKERR mField.get_finite_element_entities_by_handle("GRIFFITH_FORCE_ELEMENT",
2887 ents_to_remove);
2888 ents_to_remove = subtract(ents_to_remove, crack_front_adj_tris);
2889 CHKERR mField.remove_ents_from_finite_element("GRIFFITH_FORCE_ELEMENT",
2890 ents_to_remove, verb);
2891
2892 CHKERR mField.add_ents_to_finite_element_by_type(crack_front_adj_tris, MBTRI,
2893 "GRIFFITH_FORCE_ELEMENT");
2894 CHKERR mField.build_finite_elements("GRIFFITH_FORCE_ELEMENT",
2895 &crack_front_adj_tris, verb);
2896
2897#ifdef __ANALITICAL_TRACTION__
2898 {
2899 Range tets_level;
2900 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2901 bit, mask, MBTET, tets_level);
2902 if (proc_only) {
2903 tets_level = intersect(bitProcEnts, tets_level);
2904 }
2905
2906 Range dofs_ents;
2907 dofs_ents.merge(tets_level);
2908 for (int dd = 1; dd != 3; dd++) {
2909 CHKERR moab.get_adjacencies(tets_level, dd, false, dofs_ents,
2910 moab::Interface::UNION);
2911 }
2912 Range dofs_nodes;
2913 CHKERR moab.get_connectivity(tets_level, dofs_nodes, true);
2914 Range lower_dim_ents = unite(dofs_ents, dofs_nodes);
2915
2916 CHKERR mField.add_finite_element("ANALITICAL_METERIAL_TRACTION", MF_ZERO);
2918 "ANALITICAL_METERIAL_TRACTION", "MESH_NODE_POSITIONS");
2920 "ANALITICAL_METERIAL_TRACTION", "MESH_NODE_POSITIONS");
2922 "ANALITICAL_METERIAL_TRACTION", "SPATIAL_POSITION");
2924 "ANALITICAL_METERIAL_TRACTION", "MESH_NODE_POSITIONS");
2925 MeshsetsManager *meshset_manager_ptr;
2926 ierr = mField.getInterface(meshset_manager_ptr);
2927 CHKERRQ(ierr);
2928 if (meshset_manager_ptr->checkMeshset("ANALITICAL_TRACTION")) {
2929 const CubitMeshSets *cubit_meshset_ptr;
2930 meshset_manager_ptr->getCubitMeshsetPtr("ANALITICAL_TRACTION",
2931 &cubit_meshset_ptr);
2932 Range tris;
2933 CHKERR meshset_manager_ptr->getEntitiesByDimension(
2934 cubit_meshset_ptr->getMeshsetId(), BLOCKSET, 2, tris, true);
2935 tris = intersect(tris, lower_dim_ents);
2937 tris, MBTRI, "ANALITICAL_METERIAL_TRACTION");
2938 }
2939 CHKERR mField.build_finite_elements("ANALITICAL_METERIAL_TRACTION",
2940 &lower_dim_ents);
2941 }
2942#endif //__ANALITICAL_TRACTION__
2944}
2945
2947 const BitRefLevel bit_spatial, const BitRefLevel bit_material,
2948 const BitRefLevel mask, const bool proc_only, const bool verb) {
2950
2951 auto get_prisms_level = [&](auto bit) {
2952 Range prisms_level;
2953 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2954 bit, mask, MBPRISM, prisms_level);
2955 if (proc_only) {
2956 prisms_level = intersect(bitProcEnts, prisms_level);
2957 }
2958 return prisms_level;
2959 };
2960
2961 auto get_crack_prisms = [&](auto prisms_level) {
2962 Range crack_prisms = subtract(prisms_level, contactElements);
2963 crack_prisms = subtract(crack_prisms, mortarContactElements);
2964 return crack_prisms;
2965 };
2966
2967 auto add_both_sides_fe = [&](const std::string fe_name,
2968 const std::string lambda_field_name,
2969 const std::string spatial_field, Range ents) {
2971
2973
2974 if (!ents.empty()) {
2976 lambda_field_name);
2977 CHKERR mField.modify_finite_element_add_field_row(fe_name, spatial_field);
2978
2980 lambda_field_name);
2981 CHKERR mField.modify_finite_element_add_field_col(fe_name, spatial_field);
2982
2984 lambda_field_name);
2986 spatial_field);
2987
2988 CHKERR mField.add_ents_to_finite_element_by_type(ents, MBPRISM, fe_name);
2989 }
2990
2991 CHKERR mField.build_finite_elements(fe_name, &ents);
2992
2994 };
2995
2997 CHKERR add_both_sides_fe("CLOSE_CRACK", "LAMBDA_CLOSE_CRACK",
2998 "SPATIAL_POSITION",
2999 get_crack_prisms(get_prisms_level(bit_spatial)));
3000 }
3001 CHKERR add_both_sides_fe("BOTH_SIDES_OF_CRACK", "LAMBDA_BOTH_SIDES",
3002 "MESH_NODE_POSITIONS",
3003 get_crack_prisms(get_prisms_level(bit_material)));
3004 {
3005 Range contact_prisms =
3006 intersect(get_prisms_level(bit_material), contactElements);
3007 CHKERR add_both_sides_fe("BOTH_SIDES_OF_CONTACT",
3008 "LAMBDA_BOTH_SIDES_CONTACT", "MESH_NODE_POSITIONS",
3009 contact_prisms);
3010 }
3012}
3013
3015 const BitRefLevel mask,
3016 const bool proc_only,
3017 const bool verb) {
3019 Range tets_level;
3020 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3021 bit, mask, MBTET, tets_level);
3022 if (proc_only) {
3023 tets_level = intersect(bitProcEnts, tets_level);
3024 }
3025
3026 // Add finite elements
3027 CHKERR mField.add_finite_element("SMOOTHING", MF_ZERO, verb);
3029 "MESH_NODE_POSITIONS");
3031 "MESH_NODE_POSITIONS");
3033 "MESH_NODE_POSITIONS");
3035 "SPATIAL_POSITION");
3037 "SMOOTHING", "LAMBDA_CRACKFRONT_AREA_TANGENT");
3039 "SMOOTHING", "LAMBDA_CRACKFRONT_AREA_TANGENT");
3040
3041 // remove old ents no longer used as a this element
3042 Range ents_to_remove;
3044 ents_to_remove);
3045 ents_to_remove = subtract(ents_to_remove, tets_level);
3046 CHKERR mField.remove_ents_from_finite_element("SMOOTHING", ents_to_remove);
3047
3049 "SMOOTHING");
3050 CHKERR mField.build_finite_elements("SMOOTHING", &tets_level, verb);
3051
3053}
3054
3056 std::string fe_name, const BitRefLevel bit, const BitRefLevel mask,
3057 const std::vector<int> &ids, const bool proc_only, const int verb,
3058 const bool debug) {
3060 MeshsetsManager *meshset_manager_ptr;
3062 CHKERR mField.getInterface(meshset_manager_ptr);
3063
3064 Range level_tris;
3065 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3066 bit, BitRefLevel().set(), MBTRI, level_tris);
3067 Range surface_ents = intersect(bodySkin, level_tris);
3068
3069 if (proc_only) {
3070 // Note that some elements are on this proc, but not owned by that proc,
3071 // since
3072 // some tetrahedral share lower dimension entities
3073 Range proc_tris;
3074 CHKERR moab.get_adjacencies(bitProcEnts, 2, false, proc_tris,
3075 moab::Interface::UNION);
3076 surface_ents = intersect(surface_ents, proc_tris);
3077 }
3078
3079 CHKERR mField.add_finite_element(fe_name.c_str(), MF_ZERO, verb);
3081 "MESH_NODE_POSITIONS");
3083 "MESH_NODE_POSITIONS");
3085 "MESH_NODE_POSITIONS");
3086
3087 for (std::vector<int>::const_iterator it = ids.begin(); it != ids.end();
3088 it++) {
3089 std::string field_name =
3090 "LAMBDA_SURFACE" + boost::lexical_cast<std::string>(*it);
3092 if (dit->get()->getOwnerProc() == mField.get_comm_rank()) {
3093 EntityHandle ent = dit->get()->getEnt();
3094 Range adj_tris;
3095 CHKERR moab.get_adjacencies(&ent, 1, 2, false, adj_tris);
3096 adj_tris = intersect(adj_tris, bodySkin);
3097 surface_ents.merge(adj_tris);
3098 }
3099 }
3100 // Add finite elements
3102 field_name.c_str());
3104 field_name.c_str());
3106 field_name.c_str());
3107
3108 // remove old ents no longer used as a this element
3109 Range ents_to_remove;
3111 ents_to_remove);
3112 ents_to_remove = subtract(ents_to_remove, surface_ents);
3113 CHKERR mField.remove_ents_from_finite_element(fe_name, ents_to_remove);
3115 fe_name);
3116 }
3117
3118 CHKERR mField.build_finite_elements(fe_name, &surface_ents, verb);
3119
3121}
3122
3125 const BitRefLevel mask, const bool proc_only,
3126 const int verb, const bool debug) {
3127 auto meshset_manager_ptr = mField.getInterface<MeshsetsManager>();
3128 auto bit_ref_manager_ptr = mField.getInterface<BitRefManager>();
3130
3131 auto get_edges_block_ents = [this, meshset_manager_ptr]() {
3132 EntityHandle edge_block_set = getInterface<CPMeshCut>()->getEdgesBlockSet();
3133 Range meshset_ents;
3134 if (meshset_manager_ptr->checkMeshset(edge_block_set, BLOCKSET)) {
3135 CHKERR meshset_manager_ptr->getEntitiesByDimension(
3136 edge_block_set, BLOCKSET, 1, meshset_ents, true);
3137 }
3138 return meshset_ents;
3139 };
3140
3141 auto get_skin = [this, bit_ref_manager_ptr]() {
3142 Range tets;
3143 CHKERR bit_ref_manager_ptr->getEntitiesByTypeAndRefLevel(
3144 mapBitLevel["mesh_cut"], BitRefLevel().set(), MBTET, tets);
3145 Skinner skin(&mField.get_moab());
3146 Range skin_faces;
3147 CHKERR skin.find_skin(0, tets, false, skin_faces);
3148 return skin_faces;
3149 };
3150
3151 auto intersect_with_skin_edges = [this, get_edges_block_ents, get_skin]() {
3152 Range adj_edges;
3153 CHKERR mField.get_moab().get_adjacencies(get_skin(), 1, false, adj_edges,
3154 moab::Interface::UNION);
3155 return intersect(get_edges_block_ents(), adj_edges);
3156 };
3157 Range edges_ents = intersect_with_skin_edges();
3158
3159 Range contact_faces;
3160 contact_faces.merge(contactSlaveFaces);
3161 contact_faces.merge(contactMasterFaces);
3162
3164 boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>(
3165 new FaceOrientation(false, contact_faces, mapBitLevel["mesh_cut"]));
3166
3168 mField.get_moab(), edges_ents, get_skin());
3170 mField.get_moab(), edges_ents, get_skin(), false, contactOrientation,
3171 &mField);
3172
3173 if (debug) {
3174 Range faces = get_skin();
3176 mField.get_moab(), "edges_normals.vtk", edges_ents, &faces);
3177 }
3178
3179 auto intersect_with_bit_level_edges = [bit_ref_manager_ptr, &bit,
3180 &edges_ents]() {
3181 Range level_edges;
3182 CHKERR bit_ref_manager_ptr->getEntitiesByTypeAndRefLevel(
3183 bit, BitRefLevel().set(), MBEDGE, level_edges);
3184 return intersect(edges_ents, level_edges);
3185 };
3186 edges_ents = intersect_with_bit_level_edges();
3187
3188 auto intersect_with_proc_edges = [this, &edges_ents, proc_only]() {
3189 if (proc_only) {
3190 Range proc_edges;
3191 CHKERR mField.get_moab().get_adjacencies(
3192 bitProcEnts, 1, false, proc_edges, moab::Interface::UNION);
3193 return intersect(edges_ents, proc_edges);
3194 }
3195 return edges_ents;
3196 };
3197 edges_ents = intersect_with_proc_edges();
3198 if (proc_only)
3199 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(edges_ents,
3200 QUIET);
3201
3202 CHKERR mField.add_finite_element(fe_name.c_str(), MF_ZERO, verb);
3203
3204 auto add_field = [this, &fe_name](std::string &&field_name) {
3210 };
3211 CHKERR add_field("MESH_NODE_POSITIONS");
3212 CHKERR add_field("LAMBDA_EDGE");
3213
3214 if (debug) {
3216 mField.get_moab(),
3217 "edges_normals_on_proc_" +
3218 boost::lexical_cast<std::string>(mField.get_comm_rank()) + ".vtk",
3219 edges_ents);
3220 }
3221
3222 // remove old ents no longer used as a this element
3223 Range ents_to_remove;
3224 CHKERR mField.get_finite_element_entities_by_handle(fe_name, ents_to_remove);
3225 ents_to_remove = subtract(ents_to_remove, edges_ents);
3226 CHKERR mField.remove_ents_from_finite_element(fe_name, ents_to_remove);
3227 CHKERR mField.add_ents_to_finite_element_by_type(edges_ents, MBEDGE, fe_name);
3228
3229 CHKERR mField.build_finite_elements(fe_name, &edges_ents, verb);
3230
3232}
3233
3235 std::string fe_name, const BitRefLevel bit, const BitRefLevel mask,
3236 const bool proc_only, const int verb, const bool debug) {
3238 MeshsetsManager *meshset_manager_ptr;
3240 CHKERR mField.getInterface(meshset_manager_ptr);
3241
3242 Range level_tris;
3243 CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3244 bit, mask, MBTRI, level_tris);
3245 Range crack_surface;
3246 crack_surface = intersect(oneSideCrackFaces, level_tris);
3247 if (proc_only) {
3248 Range proc_tris;
3249 CHKERR moab.get_adjacencies(bitProcEnts, 2, false, proc_tris,
3250 moab::Interface::UNION);
3251 crack_surface = intersect(crack_surface, proc_tris);
3252 }
3253 Range crack_front_adj_tris;
3254 CHKERR mField.get_moab().get_adjacencies(
3255 crackFrontNodes, 2, false, crack_front_adj_tris, moab::Interface::UNION);
3256 crack_front_adj_tris = intersect(crack_front_adj_tris, crackFaces);
3257 crack_front_adj_tris = intersect(crack_front_adj_tris, level_tris);
3258
3259 std::string field_name =
3260 "LAMBDA_SURFACE" + boost::lexical_cast<std::string>(
3261 getInterface<CPMeshCut>()->getCrackSurfaceId());
3262
3263 // Add finite elements
3264 CHKERR mField.add_finite_element(fe_name.c_str(), MF_ZERO, verb);
3266 field_name.c_str());
3268 "MESH_NODE_POSITIONS");
3270 field_name.c_str());
3272 "MESH_NODE_POSITIONS");
3274 field_name.c_str());
3276 "MESH_NODE_POSITIONS");
3277 if (proc_only)
3278 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
3279 crack_surface, QUIET);
3280
3281 {
3282 // remove ents no longer part of this element
3283 Range ents;
3285 ents = subtract(ents, crack_surface);
3286 if (!ents.empty()) {
3288 }
3289 }
3290
3292 fe_name);
3293 CHKERR mField.build_finite_elements(fe_name, &crack_surface, verb);
3294
3295 if (debug) {
3296 std::ostringstream file;
3297 file << "out_crack_faces_" << mField.get_comm_rank() << ".vtk";
3298 EntityHandle meshset;
3299 CHKERR moab.create_meshset(MESHSET_SET, meshset);
3300 CHKERR moab.add_entities(meshset, crack_surface);
3301 CHKERR moab.write_file(file.str().c_str(), "VTK", "", &meshset, 1);
3302 CHKERR moab.delete_entities(&meshset, 1);
3303 }
3304
3305 // Add CRACKFRONT_AREA_ELEM
3306 CHKERR mField.add_finite_element("CRACKFRONT_AREA_ELEM", MF_ZERO, verb);
3307 CHKERR mField.modify_finite_element_add_field_row("CRACKFRONT_AREA_ELEM",
3308 "LAMBDA_CRACKFRONT_AREA");
3309 CHKERR mField.modify_finite_element_add_field_row("CRACKFRONT_AREA_ELEM",
3310 "MESH_NODE_POSITIONS");
3311 CHKERR mField.modify_finite_element_add_field_col("CRACKFRONT_AREA_ELEM",
3312 "LAMBDA_CRACKFRONT_AREA");
3313 CHKERR mField.modify_finite_element_add_field_col("CRACKFRONT_AREA_ELEM",
3314 "MESH_NODE_POSITIONS");
3315 CHKERR mField.modify_finite_element_add_field_data("CRACKFRONT_AREA_ELEM",
3316 "MESH_NODE_POSITIONS");
3317 CHKERR mField.modify_finite_element_add_field_data("CRACKFRONT_AREA_ELEM",
3318 "LAMBDA_CRACKFRONT_AREA");
3319 // Add CRACKFRONT_AREA_TANGENT_ELEM
3320 CHKERR mField.add_finite_element("CRACKFRONT_AREA_TANGENT_ELEM", MF_ZERO,
3321 verb);
3323 "CRACKFRONT_AREA_TANGENT_ELEM", "LAMBDA_CRACKFRONT_AREA_TANGENT");
3325 "CRACKFRONT_AREA_TANGENT_ELEM", "MESH_NODE_POSITIONS");
3327 "CRACKFRONT_AREA_TANGENT_ELEM", "LAMBDA_CRACKFRONT_AREA_TANGENT");
3329 "CRACKFRONT_AREA_TANGENT_ELEM", "MESH_NODE_POSITIONS");
3331 "CRACKFRONT_AREA_TANGENT_ELEM", "MESH_NODE_POSITIONS");
3333 "CRACKFRONT_AREA_TANGENT_ELEM", "LAMBDA_CRACKFRONT_AREA_TANGENT");
3334
3335 if (proc_only)
3336 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
3337 crack_front_adj_tris, QUIET);
3338
3339 auto remeve_ents_from_fe =
3340 [this, &crack_front_adj_tris](const std::string fe_name) {
3342 Range ents;
3344 ents = subtract(ents, crack_front_adj_tris);
3345 if (!ents.empty()) {
3347 }
3349 };
3350 CHKERR remeve_ents_from_fe("CRACKFRONT_AREA_ELEM");
3351 CHKERR remeve_ents_from_fe("CRACKFRONT_AREA_TANGENT_ELEM");
3352
3353 CHKERR mField.add_ents_to_finite_element_by_type(crack_front_adj_tris, MBTRI,
3354 "CRACKFRONT_AREA_ELEM");
3356 crack_front_adj_tris, MBTRI, "CRACKFRONT_AREA_TANGENT_ELEM");
3357
3358 if (debug) {
3359 std::ostringstream file;
3360 file << "out_crack_front_adj_tris_" << mField.get_comm_rank() << ".vtk";
3361 EntityHandle meshset;
3362 CHKERR moab.create_meshset(MESHSET_SET, meshset);
3363 CHKERR moab.add_entities(meshset, crack_front_adj_tris);
3364 CHKERR moab.write_file(file.str().c_str(), "VTK", "", &meshset, 1);
3365 CHKERR moab.delete_entities(&meshset, 1);
3366 }
3367
3368 CHKERR mField.build_finite_elements("CRACKFRONT_AREA_ELEM",
3369 &crack_front_adj_tris, verb);
3370 CHKERR mField.build_finite_elements("CRACKFRONT_AREA_TANGENT_ELEM",
3371 &crack_front_adj_tris, verb);
3372
3374}
3375
3377 SmartPetscObj<DM> &dm, const std::string prb_name,
3378 SmartPetscObj<DM> dm_elastic, SmartPetscObj<DM> dm_material,
3379 const BitRefLevel bit, const BitRefLevel mask,
3380 const std::vector<std::string> fe_list) {
3381 const MoFEM::Problem *prb_elastic_ptr;
3382 const MoFEM::Problem *prb_material_ptr;
3384 CHKERR DMMoFEMGetProblemPtr(dm_elastic, &prb_elastic_ptr);
3385 CHKERR DMMoFEMGetProblemPtr(dm_material, &prb_material_ptr);
3386 dm = createSmartDM(PETSC_COMM_WORLD, "MOFEM");
3387 CHKERR DMMoFEMCreateMoFEM(dm, &mField, prb_name.c_str(), bit, mask);
3388
3389 CHKERR DMMoFEMAddElement(dm, "ELASTIC");
3390 CHKERR DMMoFEMAddElement(dm, "ELASTIC_NOT_ALE");
3391 CHKERR DMMoFEMAddElement(dm, "SPRING");
3392 CHKERR DMMoFEMAddElement(dm, "CONTACT");
3393 CHKERR DMMoFEMAddElement(dm, "MORTAR_CONTACT");
3394 CHKERR DMMoFEMAddElement(dm, "CONTACT_POST_PROC");
3395 CHKERR DMMoFEMAddElement(dm, "SKIN"); // for postprocessing
3396 CHKERR DMMoFEMAddElement(dm, "ARC_LENGTH");
3397 CHKERR DMMoFEMAddElement(dm, "FORCE_FE");
3398 CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
3399 if (isPressureAle) {
3400 CHKERR DMMoFEMAddElement(dm, "PRESSURE_ALE");
3401 }
3402 if (areSpringsAle) {
3403 CHKERR DMMoFEMAddElement(dm, "SPRING_ALE");
3404 }
3405 if (!contactElements.empty() && !ignoreContact && !fixContactNodes) {
3406 CHKERR DMMoFEMAddElement(dm, "SIMPLE_CONTACT_ALE");
3407 CHKERR DMMoFEMAddElement(dm, "MAT_CONTACT");
3408 }
3409 CHKERR DMMoFEMAddElement(dm, "MATERIAL");
3410
3411 CHKERR DMMoFEMAddElement(dm, "GRIFFITH_FORCE_ELEMENT");
3412 for (std::vector<std::string>::const_iterator it = fe_list.begin();
3413 it != fe_list.end(); it++)
3414 CHKERR DMMoFEMAddElement(dm, it->c_str());
3415
3416 CHKERR DMMoFEMAddElement(dm, "CRACK_SURFACE");
3417 CHKERR DMMoFEMAddElement(dm, "CRACKFRONT_AREA_ELEM");
3418 CHKERR DMMoFEMAddElement(dm, "CRACKFRONT_AREA_TANGENT_ELEM");
3419#ifdef __ANALITICAL_TRACTION__
3420 CHKERR DMMoFEMAddElement(dm, "ANALITICAL_METERIAL_TRACTION");
3421#endif
3422 if (mField.check_finite_element("SMOOTHING"))
3423 CHKERR DMMoFEMAddElement(dm, "SMOOTHING");
3424
3425 if (mField.check_finite_element("EDGE"))
3426 CHKERR DMMoFEMAddElement(dm, "EDGE");
3427
3428 if (mField.check_finite_element("BOTH_SIDES_OF_CRACK"))
3429 CHKERR DMMoFEMAddElement(dm, "BOTH_SIDES_OF_CRACK");
3430
3431 if (mField.check_finite_element("BOTH_SIDES_OF_CONTACT"))
3432 CHKERR DMMoFEMAddElement(dm, "BOTH_SIDES_OF_CONTACT");
3433
3434 CHKERR DMMoFEMAddRowCompositeProblem(dm, prb_elastic_ptr->getName().c_str());
3435 CHKERR DMMoFEMAddRowCompositeProblem(dm, prb_material_ptr->getName().c_str());
3436 CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
3437 CHKERR DMSetUp(dm);
3438
3439 // Create material section. This is used by field split solver and snes
3440 // monitor
3441 PetscSection section;
3442 CHKERR mField.getInterface<ISManager>()->sectionCreate(prb_name, &section);
3443 CHKERR DMSetDefaultSection(dm, section);
3444 CHKERR DMSetDefaultGlobalSection(dm, section);
3445 // CHKERR PetscSectionView(section,PETSC_VIEWER_STDOUT_WORLD);
3446 CHKERR PetscSectionDestroy(&section);
3447
3449}
3450
3452 const std::string prb_name,
3453 const BitRefLevel bit,
3454 const BitRefLevel mask) {
3455 ProblemsManager *prb_mng_ptr;
3457 CHKERR mField.getInterface(prb_mng_ptr);
3459 dm = createSmartDM(PETSC_COMM_WORLD, "MOFEM");
3460 CHKERR DMMoFEMCreateMoFEM(dm, &mField, prb_name.c_str(), bit, mask);
3461 CHKERR DMMoFEMAddElement(dm, "ELASTIC");
3462 CHKERR DMMoFEMAddElement(dm, "SPRING");
3463 CHKERR DMMoFEMAddElement(dm, "CONTACT");
3464 CHKERR DMMoFEMAddElement(dm, "MORTAR_CONTACT");
3465 CHKERR DMMoFEMAddElement(dm, "CONTACT_POST_PROC");
3466 CHKERR DMMoFEMAddElement(dm, "SKIN");
3467 CHKERR DMMoFEMAddElement(dm, "ARC_LENGTH");
3468 CHKERR DMMoFEMAddElement(dm, "FORCE_FE");
3469 CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
3470#ifdef __ANALITICAL_TRACTION__
3471 CHKERR DMMoFEMAddElement(dm, "ANALITICAL_TRACTION");
3472#endif
3473
3474 CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
3475 mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
3476
3477 CHKERR DMSetUp(dm);
3478
3479 // Create material section. This is used by field split solver and snes
3480 // monitor
3481 PetscSection section;
3482 CHKERR mField.getInterface<ISManager>()->sectionCreate(prb_name, &section);
3483 CHKERR DMSetDefaultSection(dm, section);
3484 CHKERR DMSetDefaultGlobalSection(dm, section);
3485 // CHKERR PetscSectionView(section,PETSC_VIEWER_STDOUT_WORLD);
3486 CHKERR PetscSectionDestroy(&section);
3487
3488 mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
3489
3490 // Setting arc-length control
3491 const MoFEM::Problem *problem_ptr;
3492 CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
3493 arcCtx = boost::make_shared<ArcLengthCtx>(mField, problem_ptr->getName(),
3494 "LAMBDA_ARC_LENGTH");
3495 CHKERR arcCtx->setAlphaBeta(0, 1);
3496 auto arc_snes_ctx = boost::make_shared<ArcLengthSnesCtx>(
3497 mField, problem_ptr->getName(), arcCtx);
3498 CHKERR DMMoFEMSetSnesCtx(dm, arc_snes_ctx);
3499
3501}
3502
3504 SmartPetscObj<DM> &dm, const std::string prb_name, const BitRefLevel bit,
3505 const BitRefLevel mask) {
3506 ProblemsManager *prb_mng_ptr;
3508 CHKERR mField.getInterface(prb_mng_ptr);
3510 dm = createSmartDM(PETSC_COMM_WORLD, "MOFEM");
3511 CHKERR DMMoFEMCreateMoFEM(dm, &mField, prb_name.c_str(), bit, mask);
3512 CHKERR DMMoFEMAddElement(dm, "ELASTIC");
3513 CHKERR DMMoFEMAddElement(dm, "SPRING");
3514 CHKERR DMMoFEMAddElement(dm, "CONTACT");
3515 CHKERR DMMoFEMAddElement(dm, "MORTAR_CONTACT");
3516 CHKERR DMMoFEMAddElement(dm, "CONTACT_POST_PROC");
3517 CHKERR DMMoFEMAddElement(dm, "SKIN");
3518 CHKERR DMMoFEMAddElement(dm, "ARC_LENGTH");
3519 CHKERR DMMoFEMAddElement(dm, "FORCE_FE");
3520 CHKERR DMMoFEMAddElement(dm, "PRESSURE_FE");
3521 CHKERR DMMoFEMAddElement(dm, "CLOSE_CRACK");
3522#ifdef __ANALITICAL_TRACTION__
3523 CHKERR DMMoFEMAddElement(dm, "ANALITICAL_TRACTION");
3524#endif
3525
3526 CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
3527 mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
3528
3529 CHKERR DMSetUp(dm);
3530
3531 auto remove_higher_dofs_than_approx_dofs_from_eigen_elastic = [&]() {
3533 Range ents;
3534 CHKERR mField.get_moab().get_entities_by_handle(0, ents, true);
3535 auto prb_mng = mField.getInterface<ProblemsManager>();
3536 CHKERR prb_mng->removeDofsOnEntities("EIGEN_ELASTIC", "SPATIAL_POSITION",
3537 ents, 0, 3, approxOrder + 1, 100);
3539 };
3540 CHKERR remove_higher_dofs_than_approx_dofs_from_eigen_elastic();
3541
3542 // Create material section. This is used by field split solver and snes
3543 // monitor
3544 PetscSection section;
3545 CHKERR mField.getInterface<ISManager>()->sectionCreate(prb_name, &section);
3546 CHKERR DMSetDefaultSection(dm, section);
3547 CHKERR DMSetDefaultGlobalSection(dm, section);
3548 // CHKERR PetscSectionView(section,PETSC_VIEWER_STDOUT_WORLD);
3549 CHKERR PetscSectionDestroy(&section);
3550
3551 mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
3552
3553 // Setting arc-length control
3554 const MoFEM::Problem *problem_ptr;
3555 CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
3556 arcEigenCtx = boost::make_shared<ArcLengthCtx>(mField, problem_ptr->getName(),
3557 "LAMBDA_ARC_LENGTH");
3558 CHKERR arcEigenCtx->setAlphaBeta(0, 1);
3559 auto arc_snes_ctx = boost::make_shared<ArcLengthSnesCtx>(
3560 mField, problem_ptr->getName(), arcEigenCtx);
3561 CHKERR DMMoFEMSetSnesCtx(dm, arc_snes_ctx);
3562
3564}
3565
3567 const std::string prb_name,
3568 const BitRefLevel bit,
3569 const BitRefLevel mask) {
3570 ProblemsManager *prb_mng_ptr;
3571 MeshsetsManager *meshset_manager_ptr;
3573 CHKERR mField.getInterface(prb_mng_ptr);
3574 CHKERR mField.getInterface(meshset_manager_ptr);
3576#ifdef __ANALITICAL_DISPLACEMENT__
3577 if (meshset_manager_ptr->checkMeshset("ANALITICAL_DISP")) {
3578 const CubitMeshSets *cubit_meshset_ptr;
3579 meshset_manager_ptr->getCubitMeshsetPtr("ANALITICAL_DISP",
3580 &cubit_meshset_ptr);
3581 Range tris;
3582 CHKERR meshset_manager_ptr->getEntitiesByDimension(
3583 cubit_meshset_ptr->getMeshsetId(), BLOCKSET, 2, tris, true);
3584 if (!tris.empty()) {
3585 dm = createSmartDM(PETSC_COMM_WORLD, "MOFEM");
3586 CHKERR DMMoFEMCreateMoFEM(dm, &mField, prb_name.c_str(), bit, mask);
3587 CHKERR DMMoFEMAddElement(dm, "BC_FE");
3588 CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
3589 mField.getInterface<ProblemsManager>()->synchroniseProblemEntities =
3590 PETSC_TRUE;
3591 CHKERR DMSetUp(dm);
3592 mField.getInterface<ProblemsManager>()->synchroniseProblemEntities =
3593 PETSC_FALSE;
3594 }
3595 }
3596#endif //__ANALITICAL_DISPLACEMENT__
3598}
3599
3601 SmartPetscObj<DM> &dm, const std::string prb_name, const BitRefLevel bit,
3602 const BitRefLevel mask, const std::vector<std::string> fe_list,
3603 const bool debug) {
3606 dm = createSmartDM(PETSC_COMM_WORLD, "MOFEM");
3607 CHKERR DMMoFEMCreateMoFEM(dm, &mField, prb_name.c_str(), bit, mask);
3608 CHKERR DMMoFEMAddElement(dm, "GRIFFITH_FORCE_ELEMENT");
3609
3610 for (auto &fe : fe_list)
3611 CHKERR DMMoFEMAddElement(dm, fe.c_str());
3612
3613 CHKERR DMMoFEMAddElement(dm, "CRACK_SURFACE");
3614 CHKERR DMMoFEMAddElement(dm, "CRACKFRONT_AREA_ELEM");
3615 CHKERR DMMoFEMAddElement(dm, "CRACKFRONT_AREA_TANGENT_ELEM");
3616#ifdef __ANALITICAL_TRACTION__
3617 CHKERR DMMoFEMAddElement(dm, "ANALITICAL_METERIAL_TRACTION");
3618#endif
3619 if (mField.check_finite_element("SMOOTHING"))
3620 CHKERR DMMoFEMAddElement(dm, "SMOOTHING");
3621
3622 if (mField.check_finite_element("EDGE"))
3623 CHKERR DMMoFEMAddElement(dm, "EDGE");
3624
3625 if (mField.check_finite_element("BOTH_SIDES_OF_CRACK"))
3626 CHKERR DMMoFEMAddElement(dm, "BOTH_SIDES_OF_CRACK");
3627
3628 if (mField.check_finite_element("BOTH_SIDES_OF_CONTACT"))
3629 CHKERR DMMoFEMAddElement(dm, "BOTH_SIDES_OF_CONTACT");
3630
3631 CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
3632 CHKERR DMSetUp(dm);
3633
3634 // Create crack propagation section
3635 PetscSection section;
3636 CHKERR mField.getInterface<ISManager>()->sectionCreate(prb_name, &section);
3637 CHKERR DMSetDefaultSection(dm, section);
3638 CHKERR DMSetDefaultGlobalSection(dm, section);
3639 // CHKERR PetscSectionView(section,PETSC_VIEWER_STDOUT_WORLD);
3640 CHKERR PetscSectionDestroy(&section);
3641
3642 if (debug) {
3643 std::vector<std::string> fe_lists;
3644 const MoFEM::Problem *problem_ptr;
3645 CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
3646 const FiniteElement_multiIndex *fe_ptr;
3648 for (FiniteElement_multiIndex::iterator fit = fe_ptr->begin();
3649 fit != fe_ptr->end(); fit++) {
3650 if ((fit->get()->getId() & problem_ptr->getBitFEId()).any()) {
3651 std::string fe_name = fit->get()->getName();
3652 EntityHandle meshset;
3653 CHKERR mField.getInterface<ProblemsManager>()->getFEMeshset(
3654 prb_name, fe_name, &meshset);
3655 CHKERR mField.get_moab().write_file(
3656 (prb_name + "_" + fe_name + ".h5m").c_str(), "MOAB",
3657 "PARALLEL=WRITE_PART", &meshset, 1);
3658 CHKERR mField.get_moab().delete_entities(&meshset, 1);
3659 }
3660 }
3661 }
3663}
3664
3666 SmartPetscObj<DM> &dm, SmartPetscObj<DM> dm_material,
3667 const std::string prb_name, const int verb) {
3669 dm = createSmartDM(PETSC_COMM_WORLD, "MOFEM");
3670 CHKERR DMMoFEMSetVerbosity(dm, verb);
3671 CHKERR DMMoFEMCreateSubDM(dm, dm_material, prb_name.c_str());
3672 CHKERR DMMoFEMAddSubFieldRow(dm, "MESH_NODE_POSITIONS");
3673 CHKERR DMMoFEMAddSubFieldCol(dm, "MESH_NODE_POSITIONS");
3674 CHKERR DMMoFEMAddElement(dm, "MATERIAL");
3675 if (mField.check_finite_element("SMOOTHING")) {
3676 CHKERR DMMoFEMAddElement(dm, "SMOOTHING");
3677 }
3678 CHKERR DMMoFEMAddElement(dm, "GRIFFITH_FORCE_ELEMENT");
3679#ifdef __ANALITICAL_TRACTION__
3680 CHKERR DMMoFEMAddElement(dm, "ANALITICAL_METERIAL_TRACTION");
3681#endif
3682 CHKERR DMMoFEMSetSquareProblem(dm, PETSC_TRUE);
3683 CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
3684 CHKERR DMSetUp(dm);
3686}
3687
3689 SmartPetscObj<DM> &dm, SmartPetscObj<DM> dm_material,
3690 const std::string prb_name, const std::vector<int> surface_ids,
3691 std::vector<std::string> fe_list, const int verb) {
3693 dm = createSmartDM(PETSC_COMM_WORLD, "MOFEM");
3694 CHKERR DMMoFEMSetVerbosity(dm, verb);
3695 CHKERR DMMoFEMCreateSubDM(dm, dm_material, prb_name.c_str());
3696 CHKERR DMMoFEMSetSquareProblem(dm, PETSC_FALSE);
3697 for (auto id : surface_ids) {
3699 dm, ("LAMBDA_SURFACE" + boost::lexical_cast<std::string>(id)).c_str());
3700 }
3701 CHKERR DMMoFEMAddSubFieldCol(dm, "MESH_NODE_POSITIONS");
3702 for (auto fe : fe_list) {
3703 CHKERR DMMoFEMAddElement(dm, fe.c_str());
3704 }
3705 auto get_edges_block_set = [this]() {
3706 return getInterface<CPMeshCut>()->getEdgesBlockSet();
3707 };
3708 if (get_edges_block_set()) {
3709 CHKERR DMMoFEMAddSubFieldRow(dm, "LAMBDA_EDGE");
3710 CHKERR DMMoFEMAddElement(dm, "EDGE");
3711 }
3712 CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
3713 CHKERR DMSetUp(dm);
3715}
3716
3718 SmartPetscObj<DM> &dm, SmartPetscObj<DM> dm_material,
3719 const std::string prb_name, const bool verb) {
3721 dm = createSmartDM(PETSC_COMM_WORLD, "MOFEM");
3722 CHKERR DMMoFEMSetVerbosity(dm, verb);
3723 CHKERR DMMoFEMCreateSubDM(dm, dm_material, prb_name.c_str());
3724 CHKERR DMMoFEMAddSubFieldRow(dm, "LAMBDA_CRACKFRONT_AREA");
3725 CHKERR DMMoFEMAddSubFieldRow(dm, "LAMBDA_CRACKFRONT_AREA_TANGENT");
3726 CHKERR DMMoFEMAddSubFieldCol(dm, "MESH_NODE_POSITIONS");
3727 CHKERR DMMoFEMAddElement(dm, "CRACKFRONT_AREA_ELEM");
3728 CHKERR DMMoFEMAddElement(dm, "CRACKFRONT_AREA_TANGENT_ELEM");
3729 CHKERR DMMoFEMSetSquareProblem(dm, PETSC_FALSE);
3730 CHKERR DMMoFEMSetIsPartitioned(dm, PETSC_TRUE);
3731 CHKERR DMSetUp(dm);
3733}
3734
3736CrackPropagation::assembleElasticDM(const std::string mwls_stress_tag_name,
3737 const int verb, const bool debug) {
3739
3740 // Create elastic element data structure
3741 elasticFe = boost::make_shared<NonlinearElasticElement>(mField, ELASTIC_TAG);
3742
3743 // Create elastic element finite element instance for residuals. Note that
3744 // this element approx. singularity at crack front.
3745 feRhs = boost::make_shared<CrackFrontElement>(
3748
3749 // Create finite element instance for assembly of tangent matrix
3750 feLhs = boost::make_shared<CrackFrontElement>(
3753
3754 // Arbitrary lagrangian formulation, mesh node positions are taken into
3755 // account by operators.
3756 feRhs->meshPositionsFieldName = "NONE";
3757 feLhs->meshPositionsFieldName = "NONE";
3758 feRhs->addToRule = 0;
3759 feLhs->addToRule = 0;
3760
3761 // Create operators to calculate finite element matrices for elastic element
3762 onlyHooke = true;
3764 boost::shared_ptr<map<int, BlockData>> block_sets_ptr(
3765 elasticFe, &elasticFe->setOfBlocks);
3766 {
3768 mField, BLOCKSET | MAT_ELASTICSET, it)) {
3769
3770 // Get data from block
3771 Mat_Elastic mydata;
3772 CHKERR it->getAttributeDataStructure(mydata);
3773 int id = it->getMeshsetId();
3774 EntityHandle meshset = it->getMeshset();
3775 CHKERR mField.get_moab().get_entities_by_type(
3776 meshset, MBTET, elasticFe->setOfBlocks[id].tEts, true);
3777 elasticFe->setOfBlocks[id].iD = id;
3778 elasticFe->setOfBlocks[id].E = mydata.data.Young;
3779 elasticFe->setOfBlocks[id].PoissonRatio = mydata.data.Poisson;
3780
3781 // Print material properties of blocks after being read
3782 CHKERR PetscPrintf(PETSC_COMM_WORLD, "\nMaterial block %d \n", id);
3783 CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tYoung's modulus %6.4e \n",
3784 mydata.data.Young);
3785 CHKERR PetscPrintf(PETSC_COMM_WORLD, "\tPoisson's ratio %6.4e \n\n",
3786 mydata.data.Poisson);
3787
3788 // Set default material to elastic blocks. Create material instances
3789 switch (defaultMaterial) {
3790 case HOOKE:
3791 elasticFe->setOfBlocks[id].materialDoublePtr =
3792 boost::make_shared<Hooke<double>>();
3793 elasticFe->setOfBlocks[id].materialAdoublePtr =
3794 boost::make_shared<Hooke<adouble>>();
3795 break;
3796 case KIRCHHOFF:
3797 elasticFe->setOfBlocks[id].materialDoublePtr = boost::make_shared<
3799 double>>();
3800 elasticFe->setOfBlocks[id].materialAdoublePtr = boost::make_shared<
3802 adouble>>();
3803 onlyHooke = false;
3804 break;
3805 case NEOHOOKEAN:
3806 elasticFe->setOfBlocks[id].materialDoublePtr =
3807 boost::make_shared<NeoHookean<double>>();
3808 elasticFe->setOfBlocks[id].materialAdoublePtr =
3809 boost::make_shared<NeoHookean<adouble>>();
3810 onlyHooke = false;
3811 break;
3812 case BONEHOOKE:
3813 elasticFe->setOfBlocks[id].materialDoublePtr =
3814 boost::make_shared<Hooke<double>>();
3815 elasticFe->setOfBlocks[id].materialAdoublePtr =
3816 boost::make_shared<Hooke<adouble>>();
3817 break;
3818 default:
3819 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
3820 "Material type not yet implemented");
3821 break;
3822 }
3823 }
3824 if (onlyHookeFromOptions == PETSC_FALSE)
3825 onlyHooke = false;
3826
3827 elasticFe->commonData.spatialPositions = "SPATIAL_POSITION";
3828 elasticFe->commonData.meshPositions = "MESH_NODE_POSITIONS";
3829
3830 auto data_hooke_element_at_pts =
3831 boost::make_shared<HookeElement::DataAtIntegrationPts>();
3832
3833 // calculate material positions at integration points
3834 auto mat_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
3835 feRhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
3836 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
3837 boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
3838 if (residualStressBlock != -1 || densityMapBlock != -1) {
3839 mat_grad_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
3840 feRhs->getOpPtrVector().push_back(
3841 new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
3842 mat_grad_pos_at_pts_ptr));
3843 }
3844
3845 /// Rhs
3846 if (!onlyHooke) {
3847 // Calculate spatial positions and gradients (of deformation) at
3848 // integration pts. This operator takes into account singularity at crack
3849 // front
3850 feRhs->getOpPtrVector().push_back(new OpGetCrackFrontCommonDataAtGaussPts(
3851 "SPATIAL_POSITION", elasticFe->commonData, feRhs->singularElement,
3852 feRhs->invSJac));
3853 // Calculate material positions and gradients at integration pts.
3854 feRhs->getOpPtrVector().push_back(
3856 "MESH_NODE_POSITIONS", elasticFe->commonData));
3857 // Transfrom base functions to create singularity at crack front
3858 feRhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
3859 feRhs->singularElement, feRhs->detS, feRhs->invSJac));
3860 // Iterate over material blocks
3861 map<int, NonlinearElasticElement::BlockData>::iterator sit =
3862 elasticFe->setOfBlocks.begin();
3863 for (; sit != elasticFe->setOfBlocks.end(); sit++) {
3864 // Evaluate stress at integration pts.
3865 feRhs->getOpPtrVector().push_back(
3867 "SPATIAL_POSITION", sit->second, elasticFe->commonData,
3868 ELASTIC_TAG, false, true, false));
3869 // Assemble internal forces for right hand side
3870 feRhs->getOpPtrVector().push_back(
3872 "SPATIAL_POSITION", sit->second, elasticFe->commonData));
3873 }
3874 } else {
3875 feRhs->getOpPtrVector().push_back(
3877 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
3878 mat_grad_pos_at_pts_ptr = data_hooke_element_at_pts->HMat;
3879
3880 // Evaluate density at integration points fro material like bone
3881 if (defaultMaterial == BONEHOOKE) {
3882 if (!mwlsApprox) {
3883 // Create instance for moving least square approximation
3884 mwlsApprox = boost::make_shared<MWLSApprox>(mField);
3885 // Load mesh with stresses
3886 CHKERR mwlsApprox->loadMWLSMesh(mwlsApproxFile.c_str());
3887 MeshsetsManager *block_manager_ptr;
3888 CHKERR mField.getInterface(block_manager_ptr);
3889 CHKERR block_manager_ptr->getEntitiesByDimension(
3890 densityMapBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
3891 Tag th_rho;
3892 if (mwlsApprox->mwlsMoab.tag_get_handle(mwlsRhoTagName.c_str(),
3893 th_rho) != MB_SUCCESS) {
3894 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3895 "Density tag name %s cannot be found. Please "
3896 "provide the correct name.",
3897 mwlsRhoTagName.c_str());
3898 }
3899 CHKERR mwlsApprox->mwlsMoab.tag_get_handle(mwlsRhoTagName.c_str(),
3900 th_rho);
3901 CHKERR mwlsApprox->getValuesToNodes(th_rho);
3902 } else {
3903 mwlsApprox->tetsInBlock.clear();
3904 CHKERR mField.getInterface<MeshsetsManager>()->getEntitiesByDimension(
3905 densityMapBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
3906 }
3907 // Calculate spatial positions and gradients (of deformation) at
3908 // integration pts. This operator takes into account singularity at
3909 // crack front
3910 feRhs->getOpPtrVector().push_back(
3912 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
3913 feRhs->singularElement, feRhs->invSJac));
3914 // Evaluate density at integration pts.
3915 if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
3916 feRhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSRhoAtGaussPts(
3917 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs, mwlsApprox,
3918 mwlsRhoTagName, true, true));
3919 else {
3920 feRhs->getOpPtrVector().push_back(
3922 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs,
3923 mwlsApprox));
3924 feRhs->getOpPtrVector().push_back(
3926 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs,
3927 mwlsApprox, mwlsRhoTagName, true, false));
3928 }
3929
3930 feRhs->getOpPtrVector().push_back(
3931 new HookeElement::OpCalculateStiffnessScaledByDensityField(
3932 "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
3933 data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone,
3934 rHo0));
3935 feRhs->getOpPtrVector().push_back(
3936 new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
3937 "SPATIAL_POSITION",
3938 data_hooke_element_at_pts));
3939 feRhs->getOpPtrVector().push_back(
3940 new HookeElement::OpCalculateStress<1>("SPATIAL_POSITION",
3941 "SPATIAL_POSITION",
3942 data_hooke_element_at_pts));
3943 } else {
3944 feRhs->getOpPtrVector().push_back(
3945 new HookeElement::OpCalculateHomogeneousStiffness<0>(
3946 "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
3947 data_hooke_element_at_pts));
3948 // Calculate spatial positions and gradients (of deformation) at
3949 // integration pts. This operator takes into account singularity at
3950 // crack front
3951 feRhs->getOpPtrVector().push_back(
3953 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
3954 feRhs->singularElement, feRhs->invSJac));
3955 feRhs->getOpPtrVector().push_back(
3956 new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
3957 "SPATIAL_POSITION",
3958 data_hooke_element_at_pts));
3959 feRhs->getOpPtrVector().push_back(
3960 new HookeElement::OpCalculateStress<0>("SPATIAL_POSITION",
3961 "SPATIAL_POSITION",
3962 data_hooke_element_at_pts));
3963 }
3964
3965 feRhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
3966 feRhs->singularElement, feRhs->detS, feRhs->invSJac));
3967 feRhs->getOpPtrVector().push_back(new HookeElement::OpAleRhs_dx(
3968 "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts));
3969 }
3970
3971 // Add operators if internal stresses are present
3972 if (residualStressBlock != -1) {
3973 if (!mwlsApprox) {
3974 // Create instance for moving least square approximation
3975 mwlsApprox = boost::make_shared<MWLSApprox>(mField);
3976 // Load mesh with stresses
3977 CHKERR mwlsApprox->loadMWLSMesh(mwlsApproxFile.c_str());
3978 MeshsetsManager *block_manager_ptr;
3979 CHKERR mField.getInterface(block_manager_ptr);
3980 CHKERR block_manager_ptr->getEntitiesByDimension(
3981 residualStressBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock, true);
3982 Tag th;
3983 if (mwlsApprox->mwlsMoab.tag_get_handle(mwls_stress_tag_name.c_str(),
3984 th) != MB_SUCCESS) {
3985 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3986 "Internal stress tag name %s cannot be found. Please "
3987 "provide the correct name.",
3988 mwls_stress_tag_name.substr(4).c_str());
3989 }
3990 }
3991
3992 Tag th;
3993 if (mwlsApprox->mwlsMoab.tag_get_handle(mwls_stress_tag_name.c_str(),
3994 th) != MB_SUCCESS) {
3995 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
3996 "Internal stress tag name %s cannot be found. Please "
3997 "provide the correct name.",
3998 mwls_stress_tag_name.substr(4).c_str());
3999 }
4000 CHKERR mwlsApprox->mwlsMoab.tag_get_handle(mwls_stress_tag_name.c_str(),
4001 th);
4002 CHKERR mwlsApprox->getValuesToNodes(th);
4003
4004 // Evaluate internal stresses at integration pts.
4005 if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
4006 feRhs->getOpPtrVector().push_back(
4008 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs, mwlsApprox,
4009 mwls_stress_tag_name, false));
4010 else {
4011 feRhs->getOpPtrVector().push_back(
4013 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs,
4014 mwlsApprox));
4015 feRhs->getOpPtrVector().push_back(
4017 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feRhs, mwlsApprox,
4018 mwls_stress_tag_name, false));
4019 }
4020
4021 // Assemble internall stresses forces
4022 feRhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSSpatialStressRhs(
4023 mat_grad_pos_at_pts_ptr, mwlsApprox));
4025 feRhs->getOpPtrVector().push_back(
4026 new HookeElement::OpAnalyticalInternalAleStrain_dx<0>(
4027 "SPATIAL_POSITION", data_hooke_element_at_pts,
4029 boost::shared_ptr<MatrixDouble>(
4030 mwlsApprox, &mwlsApprox->mwlsMaterialPositions)));
4031 }
4032 }
4033
4034 /// Lhs
4035 if (!onlyHooke) {
4036 // Calculate spatial positions and gradients (of deformation) at
4037 // integration pts. This operator takes into account singularity at crack
4038 // front
4039 feLhs->getOpPtrVector().push_back(new OpGetCrackFrontCommonDataAtGaussPts(
4040 "SPATIAL_POSITION", elasticFe->commonData, feLhs->singularElement,
4041 feLhs->invSJac));
4042 // Calculate material positions and gradients at integration pts.
4043 feLhs->getOpPtrVector().push_back(
4045 "MESH_NODE_POSITIONS", elasticFe->commonData));
4046 // Transform base functions to get singular base functions at crack front
4047 feLhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
4048 feLhs->singularElement, feLhs->detS, feLhs->invSJac));
4049 // Iterate over material blocks
4050 map<int, NonlinearElasticElement::BlockData>::iterator sit =
4051 elasticFe->setOfBlocks.begin();
4052 for (; sit != elasticFe->setOfBlocks.end(); sit++) {
4053 // Evaluate stress at integration pts
4054 feLhs->getOpPtrVector().push_back(
4056 "SPATIAL_POSITION", sit->second, elasticFe->commonData,
4057 ELASTIC_TAG, true, true, false));
4058 // Assemble tanget matrix, derivative of spatial poisotions
4059 feLhs->getOpPtrVector().push_back(
4061 "SPATIAL_POSITION", "SPATIAL_POSITION", sit->second,
4062 elasticFe->commonData));
4063 }
4064 } else {
4065 // calculate material positions at integration points
4066 feLhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
4067 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
4068 boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
4069 if (residualStressBlock != -1 || densityMapBlock != -1) {
4070 mat_grad_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
4071 feLhs->getOpPtrVector().push_back(
4072 new OpCalculateVectorFieldGradient<3, 3>("MESH_NODE_POSITIONS",
4073 mat_grad_pos_at_pts_ptr));
4074 }
4075
4076 feLhs->getOpPtrVector().push_back(
4078 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
4079 mat_grad_pos_at_pts_ptr = data_hooke_element_at_pts->HMat;
4080
4081 if (defaultMaterial == BONEHOOKE) {
4082 if (!mwlsApprox) {
4083 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
4084 "mwlsApprox not allocated");
4085 } else {
4086 mwlsApprox->tetsInBlock.clear();
4087 CHKERR mField.getInterface<MeshsetsManager>()->getEntitiesByDimension(
4088 densityMapBlock, BLOCKSET, 3, mwlsApprox->tetsInBlock,
4089 true); // FIXME: do we still need this?
4090 }
4091
4092 // Calculate spatial positions and gradients (of deformation) at
4093 // integration pts. This operator takes into account singularity at
4094 // crack front
4095 feLhs->getOpPtrVector().push_back(
4097 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
4098 feLhs->singularElement, feLhs->invSJac));
4099
4100 // Evaluate density at integration pts.
4101 if (mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
4102 feLhs->getOpPtrVector().push_back(new MWLSApprox::OpMWLSRhoAtGaussPts(
4103 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feLhs, mwlsApprox,
4104 mwlsRhoTagName, true, true));
4105 else {
4106 feLhs->getOpPtrVector().push_back(
4108 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feLhs,
4109 mwlsApprox));
4110 feLhs->getOpPtrVector().push_back(
4112 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, feLhs,
4113 mwlsApprox, mwlsRhoTagName, true, false));
4114 }
4115
4116 feLhs->getOpPtrVector().push_back(
4117 new HookeElement::OpCalculateStiffnessScaledByDensityField(
4118 "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
4119 data_hooke_element_at_pts, mwlsApprox->rhoAtGaussPts, nBone,
4120 rHo0));
4121 feLhs->getOpPtrVector().push_back(
4122 new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
4123 "SPATIAL_POSITION",
4124 data_hooke_element_at_pts));
4125 feLhs->getOpPtrVector().push_back(
4126 new HookeElement::OpCalculateStress<1>("SPATIAL_POSITION",
4127 "SPATIAL_POSITION",
4128 data_hooke_element_at_pts));
4129
4130 feLhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
4131 feLhs->singularElement, feLhs->detS, feLhs->invSJac));
4132 feLhs->getOpPtrVector().push_back(new HookeElement::OpAleLhs_dx_dx<1>(
4133 "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts));
4134
4135 } else {
4136 feLhs->getOpPtrVector().push_back(
4137 new HookeElement::OpCalculateHomogeneousStiffness<0>(
4138 "SPATIAL_POSITION", "SPATIAL_POSITION", block_sets_ptr,
4139 data_hooke_element_at_pts));
4140 // Calculate spatial positions and gradients (of deformation) at
4141 // integration pts. This operator takes into account singularity at
4142 // crack front
4143 feLhs->getOpPtrVector().push_back(
4145 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
4146 feLhs->singularElement, feLhs->invSJac));
4147 feLhs->getOpPtrVector().push_back(
4148 new HookeElement::OpCalculateStrainAle("SPATIAL_POSITION",
4149 "SPATIAL_POSITION",
4150 data_hooke_element_at_pts));
4151 feLhs->getOpPtrVector().push_back(
4152 new HookeElement::OpCalculateStress<0>("SPATIAL_POSITION",
4153 "SPATIAL_POSITION",
4154 data_hooke_element_at_pts));
4155 // Transform base functions to get singular base functions at crack
4156 // front
4157 feLhs->getOpPtrVector().push_back(new OpTransfromSingularBaseFunctions(
4158 feLhs->singularElement, feLhs->detS, feLhs->invSJac));
4159 feLhs->getOpPtrVector().push_back(new HookeElement::OpAleLhs_dx_dx<0>(
4160 "SPATIAL_POSITION", "SPATIAL_POSITION", data_hooke_element_at_pts));
4161 }
4162 }
4163 }
4164
4165 // Assembly forces
4166
4168 boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
4169 // Surface forces
4170 {
4171 string fe_name_str = "FORCE_FE";
4172 surfaceForces->insert(fe_name_str, new NeumannForcesSurface(mField));
4174 it)) {
4175 CHKERR surfaceForces->at(fe_name_str)
4176 .addForce("SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(), true);
4177 }
4178
4179 const string block_set_force_name("FORCE");
4180 // search for block named FORCE and add its attributes to FORCE_FE element
4182 if (it->getName().compare(0, block_set_force_name.length(),
4183 block_set_force_name) == 0) {
4184 CHKERR surfaceForces->at(fe_name_str)
4185 .addForce("SPATIAL_POSITION", PETSC_NULL, (it->getMeshsetId()),
4186 true, true);
4187 }
4188 }
4189 }
4190 // assemble surface pressure
4192 boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
4193 {
4194 string fe_name_str = "PRESSURE_FE";
4195 surfacePressure->insert(fe_name_str, new NeumannForcesSurface(mField));
4196 // iterate over sidestep where pressure is applied