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