29 #include <moab/SpatialLocator.hpp>
30 #include <moab/ElemEvaluator.hpp>
33 using namespace MoFEM;
54 #include <Smoother.hpp>
55 #include <VolumeLengthQuality.hpp>
56 #include <SurfaceSlidingConstrains.hpp>
63 #include <petsc/private/pcimpl.h>
72 PetscViewerAndFormat *vf) {
73 PetscViewer viewer = vf->viewer;
78 PetscReal *lnorms, *norms;
79 PetscInt numFields,
f, pStart, pEnd, p;
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;
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]));
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,
110 for (
f = 0;
f < numFields; ++
f) {
114 (
double)PetscSqrtReal(norms[
f]));
116 CHKERR PetscViewerASCIISubtractTab(viewer, ((PetscObject)snes)->tablevel);
117 CHKERR PetscViewerPopFormat(viewer);
118 CHKERR PetscFree2(lnorms, norms);
124 namespace po = boost::program_options;
129 boost::shared_ptr<WrapMPIComm> moab_comm_wrap) {
131 std::vector<ParallelComm *> list_pcomms;
132 ParallelComm::get_all_pcomm(&moab, list_pcomms);
140 for (
auto p : list_pcomms) {
145 new ParallelComm(&moab, moab_comm_wrap->get_comm(), &pcomm_id);
151 const int from_proc,
Range &entities,
152 const bool adjacencies,
const bool tags) {
153 ErrorCode result = MB_SUCCESS;
157 auto add_verts = [&](
Range &sent_ents) {
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;
167 tmp_result = moab_tmp.get_entities_by_type(*rit, MBVERTEX, sent_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");
181 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab_tmp,
MYPCOMM_INDEX);
182 auto &procConfig = pcomm->proc_config();
183 const int MAX_BCAST_SIZE = (1 << 28);
185 ParallelComm::Buffer buff(ParallelComm::INITIAL_BUFF_SIZE);
186 buff.reset_ptr(
sizeof(
int));
187 if ((
int)procConfig.proc_rank() == from_proc) {
189 result = add_verts(entities);
192 buff.reset_ptr(
sizeof(
int));
193 result = pcomm->pack_buffer(entities, adjacencies, tags,
false, -1, &buff);
195 "Failed to compute buffer size in broadcast_entities");
196 buff.set_stored_size();
197 buff_size = buff.buff_ptr - buff.mem_ptr;
201 MPI_Bcast(&buff_size, 1, MPI_INT, from_proc, procConfig.proc_comm());
202 if (MPI_SUCCESS != success) {
209 if ((
int)procConfig.proc_rank() != from_proc)
210 buff.reserve(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) {
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);
235 result = pcomm->unpack_buffer(buff.buff_ptr,
false, from_proc, -1, dum1a,
236 dum1b, dum1p, dum2, dum2, dum3, dum4);
238 CHK_MOAB_THROW(result,
"Failed to unpack buffer in broadcast_entities");
239 std::copy(dum4.begin(), dum4.end(), range_inserter(entities));
249 CHKERRABORT(PETSC_COMM_WORLD,
ierr);
255 PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Fixing nodes options",
"none");
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",
264 ierr = PetscOptionsEnd();
270 const bool fix_small_g =
true,
271 const bool debug =
false)
const {
276 Range tets_level, prisms_level;
285 tets_level.merge(prisms_level);
288 CHKERR skin.find_skin(0, tets_level,
false, mesh_skin);
290 mesh_skin = mesh_skin.subset_by_type(MBTRI);
292 Range faces_to_remove, contact_faces;
295 faces_to_remove = subtract(faces_to_remove, contact_faces);
297 mesh_skin = subtract(mesh_skin, faces_to_remove);
299 CHKERR m_field.
get_moab().get_connectivity(mesh_skin, fix_material_nodes);
302 Range fixed_edges_nodes;
303 CHKERR m_field.
get_moab().get_connectivity(fixed_edges, fixed_edges_nodes,
306 Range crack_faces_and_fixed_edges_nodes;
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);
312 Range crack_faces_edges;
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;
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);
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(
331 CHKERR m_field.
get_moab().write_file(
"fixed_faces.vtk",
"VTK",
"",
340 auto bc_fix_nodes = [&m_field,
344 CHKERR m_field.
get_moab().get_entities_by_handle(meshset, ents,
true);
347 fix_material_nodes.merge(nodes);
352 CHKERR bc_fix_nodes(it->getMeshset());
358 CHKERR bc_fix_nodes(it->getMeshset());
366 CHKERR bc_fix_nodes(it->getMeshset());
375 fix_material_nodes.merge(contact_nodes);
381 if (
bit->getName().compare(0, 9,
"SPRING_BC") == 0) {
391 fix_material_nodes.merge(contact_nodes);
396 Range interface_nodes;
399 fix_material_nodes.merge(interface_nodes);
402 auto add_nodes_on_opposite_side_of_prism = [&m_field, &fix_material_nodes,
403 this](
bool only_contact_prisms =
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) {
413 for (
auto p : adj_prisms) {
416 CHKERR m_field.
get_moab().get_connectivity(p, conn, num_nodes,
false);
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]);
422 if (fix_material_nodes.find(conn[3 +
n]) !=
423 fix_material_nodes.end()) {
424 add_nodes.insert(conn[
n]);
427 fix_material_nodes.merge(add_nodes);
432 CHKERR add_nodes_on_opposite_side_of_prism();
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;
448 const std::map<EntityHandle, double> map_j = copy_map(cP.
mapJ);
449 const std::map<EntityHandle, double> map_g1 = copy_map(cP.
mapG1);
451 auto find_max = [](
const auto &map) {
454 max = fmax(
m.second, max);
458 const double max_j = find_max(map_j);
459 const double max_g1 = find_max(map_g1);
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);
472 "No points to move on crack front");
488 CHKERR add_nodes_on_opposite_side_of_prism(
true);
498 CrackPropagation::query_interface(boost::typeindex::type_index type_index,
503 if (type_index == boost::typeindex::type_id<CrackPropagation>()) {
508 if (type_index == boost::typeindex::type_id<CPSolvers>()) {
509 CHKERR cpSolversPtr->query_interface(type_index, iface);
513 if (type_index == boost::typeindex::type_id<CPMeshCut>()) {
514 CHKERR cpMeshCutPtr->query_interface(type_index, iface);
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)),
535 if (!LogManager::checkIfChannelExist(
"CPWorld")) {
536 auto core_log = logging::core::get();
539 LogManager::createSink(LogManager::getStrmWorld(),
"CPWorld"));
541 LogManager::createSink(LogManager::getStrmSync(),
"CPSync"));
543 LogManager::createSink(LogManager::getStrmSelf(),
"CPSelf"));
545 LogManager::setLog(
"CPWorld");
546 LogManager::setLog(
"CPSync");
547 LogManager::setLog(
"CPSelf");
554 MOFEM_LOG(
"CPWorld", Sev::noisy) <<
"CPSolve created";
556 #ifdef GIT_UM_SHA1_NAME
557 MOFEM_LOG_C(
"CPWorld", Sev::inform,
"User module git commit id %s",
564 <<
"Fracture module version " << version.
strVersion();
566 #ifdef GIT_FM_SHA1_NAME
567 MOFEM_LOG_C(
"CPWorld", Sev::inform,
"Fracture module git commit id %s",
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);
583 moabCommWorld = boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD,
false);
593 CHKERR PetscOptionsBool(
"-my_propagate_crack",
598 CHKERR PetscOptionsInt(
"-my_geom_order",
"approximation geometry order",
"",
600 CHKERR PetscOptionsScalar(
"-my_gc",
"release energy",
"",
gC, &
gC,
602 CHKERR PetscOptionsScalar(
"-beta_gc",
"heterogeneous gc beta coefficient",
604 CHKERR PetscOptionsBool(
"-my_is_partitioned",
"true if mesh is partitioned",
606 CHKERR PetscOptionsInt(
"-my_ref",
"crack tip mesh refinement level",
"",
608 CHKERR PetscOptionsInt(
"-my_ref_order",
609 "crack tip refinement approximation order level",
"",
612 MOFEM_LOG_C(
"CPWorld", Sev::inform,
"### Input parameter: -my_order %d",
616 MOFEM_LOG_C(
"CPWorld", Sev::inform,
"### Input parameter: -my_gc %6.4e",
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",
628 "if false force use of nonlinear element for hooke material",
"",
632 CHKERR PetscOptionsBool(
"-my_add_pressure_ale",
633 "if set surface pressure is considered in ALE",
"",
637 CHKERR PetscOptionsBool(
"-my_add_springs_ale",
638 "if set surface springs is considered in ALE",
"",
642 CHKERR PetscOptionsBool(
"-my_add_surface_force_ale",
643 "if set surface force is considered in ALE",
"",
648 "-my_ignore_material_surface_force",
649 "if set material forces arising from surface force are ignored",
"",
653 CHKERR PetscOptionsBool(
"-my_add_singularity",
654 "if set singularity added to crack front",
"",
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) {
668 "### Input parameter: -my_mwls_approx_file %s", file_name);
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);
675 "### Input parameter: -my_internal_stress_name %s", tag_name);
677 CHKERR PetscOptionsString(
"-my_eigen_stress_name",
678 "name of eigen stress tag",
"", tag_name,
679 tag_name, 255, &flg);
681 "### Input parameter: -my_eigen_stress_name %s", tag_name);
685 "-my_residual_stress_block",
686 "block in mechanical (i.e. cracked) mesh to which residual stress "
687 "and density mapping are applied",
690 "### Input parameter: -my_residual_stress_block %d",
693 CHKERR PetscOptionsInt(
"-my_density_block",
694 "block to which density is mapped",
"",
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);
704 "### Input parameter: -my_density_tag_name %s",
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);
718 "File configuring meshsets ( %s ) cannot be open\n",
730 MOFEM_LOG_C(
"CPWorld", Sev::inform,
"### Input parameter: -material %s",
735 CHKERR PetscOptionsScalar(
"-my_rho0",
"release energy",
"",
rHo0, &
rHo0,
737 MOFEM_LOG_C(
"CPWorld", Sev::inform,
"### Input parameter: -my_rho0 %6.4e",
740 CHKERR PetscOptionsScalar(
"-my_n_bone",
"release energy",
"",
nBone,
743 "### Input parameter: -my_n_bone %6.4e",
nBone);
747 CHKERR PetscOptionsInt(
"-nb_load_steps",
"number of load steps",
"",
750 CHKERR PetscOptionsInt(
"-nb_cut_steps",
"number of cut mesh steps",
"",
756 "### Input parameter: -load_scale %6.4e",
loadScale);
761 "-other_side_constrain",
762 "Set surface constrain on other side of crack surface",
"",
769 CHKERR PetscOptionsReal(
"-smoother_alpha",
"Control mesh smoother",
"",
771 CHKERR PetscOptionsReal(
"-my_gamma",
"Controls mesh smoother",
"",
775 "### Input parameter: -smoother_alpha %6.4e",
smootherAlpha);
782 CHKERR PetscOptionsReal(
"-arc_alpha",
"Arc length alpha",
"",
arcAlpha,
784 CHKERR PetscOptionsReal(
"-arc_beta",
"Arc length beta",
"",
arcBeta,
786 CHKERR PetscOptionsReal(
"-arc_s",
"Arc length step size",
"",
arcS, &
arcS,
789 "### Input parameter: -arc_alpha %6.4e",
arcAlpha);
791 "### Input parameter: -arc_beta %6.4e",
arcBeta);
792 MOFEM_LOG_C(
"CPWorld", Sev::inform,
"### Input parameter: -arc_s %6.4e",
805 CHKERR PetscOptionsBool(
"-my_ignore_contact",
"If true ignore contact",
807 CHKERR PetscOptionsBool(
"-my_fix_contact_nodes",
810 CHKERR PetscOptionsBool(
"-my_print_contact_state",
811 "If true print contact state",
"",
814 "-my_contact_output_integ_pts",
815 "If true output data at contact integration points",
"",
817 CHKERR PetscOptionsReal(
"-my_r_value",
"Contact regularisation parameter",
819 CHKERR PetscOptionsReal(
"-my_cn_value",
"Contact augmentation parameter",
821 CHKERR PetscOptionsInt(
"-my_contact_order",
"contact approximation order",
824 "-my_contact_lambda_order",
"contact Lagrange multipliers order",
"",
827 "### Input parameter: -my_cn_value %6.4e",
cnValue);
829 "### Input parameter: -my_contact_order %d",
contactOrder);
831 "### Input parameter: -my_cn_value %6.4e",
cnValue);
833 "### Input parameter: -my_contact_order %d",
contactOrder);
839 CHKERR PetscOptionsBool(
"-cut_mesh",
"If true mesh is cut by surface",
"",
842 CHKERR PetscOptionsReal(
"-cut_factor",
"Crack acceleration factor",
"",
847 "-run_elastic_without_crack",
"If true mesh is cut by surface",
"",
850 MOFEM_LOG_C(
"CPWorld", Sev::inform,
"### Input parameter: -cut_mesh %d",
853 "### Input parameter: -cut_factor %6.4e",
856 "### Input parameter: -run_elastic_without_crack %d",
859 CHKERR PetscOptionsInt(
"-post_proc_level",
"level of output files",
"",
862 "### Input parameter: -post_proc_level %6.4e",
postProcLevel);
866 "-add_analytical_internal_stress_operators",
867 "If true add analytical internal stress operators for mwls test",
"",
872 CHKERR PetscOptionsBool(
"-solve_eigen_problem",
"If true solve eigen problem",
877 "-use_eigen_pos_simple_contact",
878 "If true use eigen positions for matching meshes contact",
"",
882 "### Input parameter: -solve_eigen_problem %d",
885 "### Input parameter: -use_eigen_pos_simple_contact %d",
890 CHKERR PetscOptionsScalar(
"-part_weight_power",
"Partitioning weight power",
896 "-calc_mwls_coeffs_every_propagation_step",
897 "If true recaulculate MWLS coefficients every propagation step",
"",
901 ierr = PetscOptionsEnd();
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);
923 CHKERR getInterface<CPMeshCut>()->getOptions();
924 CHKERR getInterface<CPSolvers>()->getOptions();
933 CIRCULAR_PLATE_IN_INFINITE_BODY,
934 ANGLE_ANALYTICAL_MODE_1,
935 MWLS_INTERNAL_STRESS,
936 MWLS_INTERNAL_STRESS_ANALYTICAL,
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",
947 "notch_with_density"};
949 int choice_value = ANALYTICAL_MODE_1;
954 &choice_value, &flg);
955 CHKERR PetscOptionsScalar(
"-test_tol",
"test tolerance",
"", test_tol,
956 &test_tol, PETSC_NULL);
957 ierr = PetscOptionsEnd();
960 if (flg == PETSC_TRUE) {
964 switch (choice_value) {
965 case ANALYTICAL_MODE_1: {
967 for (map<EntityHandle, double>::iterator mit =
mapG1.begin();
968 mit !=
mapG1.end(); mit++) {
969 ave_g1 += mit->second;
971 ave_g1 /=
mapG1.size();
973 const double gc = 1e-5;
975 const double error = fabs(ave_g1 - gc) / gc;
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) {
982 "ANALYTICAL_MODE_1 test failed");
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;
995 "g = %6.4e gc = %6.4e tolerance = %6.4e error = %6.4e\n",
g, gc,
998 if (max_error > test_tol || max_error != max_error)
1000 "CIRCULAR_PLATE_IN_INFINITE_BODY test failed");
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;
1013 "g = %6.4e gc = %6.4e tolerance = %6.4e error = %6.4e\n",
g, gc,
1016 if (max_error > test_tol || max_error != max_error) {
1018 "ANGLE_ANALYTICAL_MODE_1 test failed %6.4e > %6.4e",
1019 max_error, test_tol);
1022 case MWLS_INTERNAL_STRESS: {
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;
1033 "g = %6.4e gc = %6.4e tolerance = %6.4e error = %6.4e\n",
g, gc,
1036 if (max_error > test_tol || max_error != max_error) {
1038 "MWLS_INTERNAL_STRESS test failed");
1041 case MWLS_INTERNAL_STRESS_ANALYTICAL: {
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;
1052 "g = %6.4e gc = %6.4e tolerance = %6.4e error = %6.4e\n",
g, gc,
1055 if (max_error > test_tol || max_error != max_error) {
1057 "MWLS_INTERNAL_STRESS_ANALYTICAL test failed");
1060 case CUT_CIRCLE_PLATE: {
1064 const double g =
m.second;
1065 max_g = fmax(
g, max_g);
1069 double error = fabs(max_g - gc) / gc;
1070 if (error > test_tol || error != error) {
1072 "CUT_CIRCLE_PLATE test failed");
1075 case NOTCH_WITH_DENSITY: {
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);
1083 const double gc = 1.6642e-04;
1085 double error = fabs(max_g - gc) / gc;
1086 if (error > test_tol || error != error) {
1088 "NOTCH_WITH_DENSITY test failed");
1099 const Range &ents) {
1101 std::ostringstream file;
1107 &meshset_out, 1, NULL, 0);
1121 map<int, Range> layers;
1124 CHKERR moab.get_adjacencies(ref_nodes, 3,
false, tets,
1125 moab::Interface::UNION);
1131 CHKERR moab.get_connectivity(tets, ref_nodes,
false);
1132 CHKERR moab.get_adjacencies(ref_nodes, 3,
false, tets,
1133 moab::Interface::UNION);
1137 layers[ll + 1] = subtract(tets, layers[ll]);
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;
1149 rval = moab.tag_get_handle(
"TETS_WEIGHT", th_tet_weight);
1150 if (
rval == MB_SUCCESS) {
1151 CHKERR moab.tag_delete(th_tet_weight);
1154 CHKERR moab.tag_get_handle(
"TETS_WEIGHT", 1, MB_TYPE_INTEGER, th_tet_weight,
1155 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
1157 for (
unsigned int ll = 0; ll != layers.size(); ll++) {
1159 CHKERR moab.tag_clear_data(th_tet_weight, layers[ll], &weight);
1172 CHKERR moab.tag_delete(th_tet_weight);
1178 Range &proc_ents,
const int verb,
1184 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
1187 Range all_proc_ents;
1190 Tag part_tag = pcomm->part_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();
1197 CHKERR moab.tag_get_data(part_tag, &*mit, 1, &part);
1199 CHKERR moab.get_entities_by_dimension(*mit, 3, proc_ents,
true);
1200 CHKERR moab.get_entities_by_handle(*mit, all_proc_ents,
true);
1203 proc_ents = intersect(proc_ents, tets);
1207 CHKERR skin.find_skin(0, tets,
false, tets_skin);
1210 Range proc_ents_skin[4];
1217 bit2_again,
BitRefLevel().set(), MBPRISM, prisms_level);
1221 Range contact_prisms_to_remove =
1223 prisms_level = subtract(prisms_level, contact_prisms_to_remove);
1225 Range face_on_prism;
1226 CHKERR moab.get_adjacencies(contact_prisms, 2,
false, face_on_prism,
1227 moab::Interface::UNION);
1229 Range tris_on_prism = face_on_prism.subset_by_type(MBTRI);
1231 Range adj_vols_to_prisms;
1232 CHKERR moab.get_adjacencies(tris_on_prism, 3,
false, adj_vols_to_prisms,
1233 moab::Interface::UNION);
1237 bit2_again,
BitRefLevel().set(), MBTET, tet_level);
1239 Range contact_tets_from_vols = adj_vols_to_prisms.subset_by_type(MBTET);
1241 contact_tets = intersect(tet_level, contact_tets_from_vols);
1244 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
1245 Range contact_tets_container = contact_tets;
1274 proc_ents.merge(contact_tets);
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);
1290 proc_ents_skin[3] = proc_ents;
1291 CHKERR skin.find_skin(0, proc_ents,
false, proc_ents_skin[2]);
1293 proc_ents_skin[2] = subtract(proc_ents_skin[2], tets_skin);
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);
1305 Range crack_faces_nodes;
1307 proc_ents_skin[0].merge(crack_faces_nodes);
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());
1321 proc_ents_skin[0].merge(arc_length_vertices);
1323 const double coords[] = {0, 0, 0};
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]);
1341 all_ents_no_on_part = subtract(
1342 all_ents_no_on_part, all_ents_no_on_part.subset_by_type(MBENTITYSET));
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());
1354 CHKERR pcomm->resolve_shared_ents(0, proc_ents, 3, -1, proc_ents_skin);
1359 std::ostringstream file_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]);
1367 CHKERR moab.write_file(file_skin.str().c_str(),
"VTK",
"", &meshset_skin,
1369 std::ostringstream file_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,
1376 CHKERR moab.delete_entities(&meshset_owned, 1);
1382 CHKERR pcomm->get_shared_entities(-1, shared_ents);
1383 std::ostringstream file_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);
1401 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
1409 std::ostringstream file_owned;
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);
1419 Tag part_tag = pcomm->part_tag();
1422 0, MBENTITYSET, &part_tag, NULL, 1, tagged_sets,
1423 moab::Interface::UNION);
1432 tris_level = intersect(tris_level, prisms_sides);
1437 for (Range::iterator mit = tagged_sets.begin(); mit != tagged_sets.end();
1440 CHKERR moab.get_entities_by_type(*mit, MBTRI, part_tris);
1441 part_tris = intersect(part_tris, tris_level);
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);
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);
1462 const bool build_fields,
const int verb,
const bool debug) {
1466 auto write_file = [&](
const std::string file_name,
Range ents) {
1479 Range tets_level, tets_level_all;
1481 bit, mask, MBTET, tets_level_all);
1483 tets_level = intersect(
bitProcEnts, tets_level_all);
1488 bit, mask, MBPRISM, prisms_level);
1490 prisms_level = intersect(prisms_level,
bitProcEnts);
1492 Range prisms_level_tris;
1494 prisms_level, 2,
false, prisms_level_tris, moab::Interface::UNION);
1496 if (
rval != MB_SUCCESS &&
rval != MB_MULTIPLE_ENTITIES_FOUND) {
1501 prisms_level_tris = intersect(prisms_level_tris,
crackFaces);
1502 Range proc_ents = tets_level;
1503 proc_ents.merge(prisms_level_tris);
1506 MB_TAG_SPARSE,
MF_ZERO, verb);
1513 MB_TAG_SPARSE,
MF_ZERO, verb);
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);
1524 CHKERR moab.get_connectivity(proc_ents, dofs_nodes,
true);
1530 ents_to_remove_from_field;
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);
1537 ents_to_remove_from_field, verb);
1539 ents_to_remove_from_field, verb);
1542 ents_to_remove_from_field, verb);
1545 auto add_spatial_field_ents = [&](
const auto field_name) {
1559 CHKERR add_spatial_field_ents(
"SPATIAL_POSITION");
1561 CHKERR add_spatial_field_ents(
"EIGEN_SPATIAL_POSITIONS");
1564 Range current_field_ents;
1566 current_field_ents);
1567 dofs_nodes = current_field_ents.subset_by_type(MBVERTEX);
1568 dofs_ents = subtract(current_field_ents, dofs_nodes);
1570 "MESH_NODE_POSITIONS", verb);
1572 auto set_spatial_field_order = [&](
const auto field_name) {
1578 CHKERR set_spatial_field_order(
"SPATIAL_POSITION");
1580 CHKERR set_spatial_field_order(
"EIGEN_SPATIAL_POSITIONS");
1597 non_ref_ents.merge(dofs_ents);
1599 Range contact_faces;
1609 bit, mask, MBTET, proc_ents);
1610 map<int, Range> layers;
1614 CHKERR moab.get_adjacencies(ref_nodes, 3,
false, tets,
1615 moab::Interface::UNION);
1616 tets = intersect(tets_level_all, 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);
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);
1635 Range contact_elements;
1639 contact_elements = intersect(contact_prisms, contact_elements);
1640 if (contact_elements.empty()) {
1642 "Expecting adjacent contact prism(s), but none found");
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) {
1650 "Expecting 2 or more adjacent contact tris, but %d found",
1651 contact_tris.size());
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);
1666 for (map<int, Range>::reverse_iterator mit = layers.rbegin();
1667 mit != layers.rend(); mit++, poo++) {
1668 mit->second = intersect(mit->second, dofs_ents);
1671 non_ref_ents = subtract(non_ref_ents, mit->second);
1683 Range cont_ents, cont_adj_ent;
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);
1727 const bool build_fields,
1736 if (vit == refined_ents_ptr->get<
Ent_mi_tag>().end()) {
1738 "Vertex of arc-length field not found");
1740 *(
const_cast<RefEntity *
>(vit->get())->getBitRefLevelPtr()) |=
bit;
1742 FieldEntity::getLocalUniqueIdCalculate(
1746 "Arc length field entity not found");
1748 MOFEM_LOG(
"CPWorld", Sev::noisy) << **fit;
1749 MOFEM_LOG_C(
"CPWorld", Sev::inform,
"Arc-Length field lambda = %3.4e",
1750 (*fit)->getEntFieldData()[0]);
1755 MB_TAG_SPARSE,
MF_ZERO, verb);
1759 std::pair<RefEntity_multiIndex::iterator, bool> p_ent =
1761 ->insert(boost::make_shared<RefEntity>(
1763 *(
const_cast<RefEntity *
>(p_ent.first->get())->getBitRefLevelPtr()) |=
1772 const bool proc_only,
1773 const bool build_fields,
1785 Range body_skin_proc_tris;
1786 body_skin_proc_tris = intersect(
bodySkin, level_tris);
1791 moab::Interface::UNION);
1792 body_skin_proc_tris = intersect(body_skin_proc_tris, proc_tris);
1809 "LAMBDA_SURFACE" + boost::lexical_cast<std::string>(
1810 getInterface<CPMeshCut>()->getSkinOfTheBodyId());
1812 CHKERR moab.get_connectivity(body_skin_proc_tris, dofs_nodes,
true);
1818 MB_TAG_SPARSE,
MF_ZERO, verb);
1820 Range ents_to_remove;
1822 ents_to_remove = subtract(ents_to_remove, dofs_nodes);
1837 const bool proc_only,
1838 const bool build_fields,
1847 CHKERR bit_ref_manager_ptr->getEntitiesByTypeAndRefLevel(
1850 if (meshset_manager_ptr->checkMeshset(
1851 getInterface<CPMeshCut>()->getEdgesBlockSet(),
BLOCKSET)) {
1853 CHKERR meshset_manager_ptr->getEntitiesByDimension(
1854 getInterface<CPMeshCut>()->getEdgesBlockSet(),
BLOCKSET, 1,
1855 meshset_ents,
true);
1856 edges_ents = intersect(meshset_ents, level_edges);
1862 moab::Interface::UNION);
1863 edges_ents = intersect(edges_ents, proc_edges);
1867 Range other_side_crack_faces_nodes;
1869 other_side_crack_faces_nodes,
true);
1870 other_side_crack_faces_nodes =
1874 CHKERR moab.get_connectivity(edges_ents, dofs_nodes,
true);
1875 dofs_nodes = subtract(dofs_nodes, other_side_crack_faces_nodes);
1882 MB_TAG_SPARSE,
MF_ZERO, verb);
1884 Range ents_to_remove;
1886 ents_to_remove = subtract(ents_to_remove, dofs_nodes);
1893 auto get_fields_multi_index = [
this]() {
1898 for (
auto field : get_fields_multi_index()) {
1899 if (field->getName().compare(0, 14,
"LAMBDA_SURFACE") == 0) {
1912 const BitRefLevel bit,
const bool proc_only,
const bool build_fields,
1913 const int verb,
const bool debug) {
1923 crack_surface_ents = intersect(crack_surface_ents, level_tris);
1928 moab::Interface::UNION);
1929 crack_surface_ents = intersect(crack_surface_ents, proc_tris);
1946 "LAMBDA_SURFACE" + boost::lexical_cast<std::string>(
1947 getInterface<CPMeshCut>()->getCrackSurfaceId());
1949 CHKERR moab.get_connectivity(crack_surface_ents, dofs_nodes,
true);
1956 MB_TAG_SPARSE,
MF_ZERO, verb);
1957 Range ents_to_remove;
1959 ents_to_remove = subtract(ents_to_remove, dofs_nodes);
1975 const bool proc_only,
const bool build_fields,
const int verb,
1979 auto add_both_sides_field = [&](
const std::string lambda_field_name,
1989 std::ostringstream file;
1990 file << lambda_field_name <<
".vtk";
2000 Range ents_to_remove;
2003 ents_to_remove = subtract(ents_to_remove, master_nodes);
2016 lambda_field_name, 1);
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) {
2035 if (
field_name.compare(0, 11,
"LAMBDA_EDGE") == 0) {
2042 auto write_file = [&](
const std::string file_name,
Range ents) {
2055 auto get_prisms_level_nodes = [&](
auto bit) {
2060 prisms_level = intersect(
bitProcEnts, prisms_level);
2062 Range prisms_level_nodes;
2065 return prisms_level_nodes;
2068 auto get_other_side_crack_faces_nodes = [&](
auto prisms_level_nodes) {
2069 Range other_side_crack_faces_nodes;
2072 other_side_crack_faces_nodes =
2073 intersect(other_side_crack_faces_nodes, prisms_level_nodes);
2074 other_side_crack_faces_nodes =
2077 return other_side_crack_faces_nodes;
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);
2084 CHKERR add_both_sides_field(
"LAMBDA_BOTH_SIDES", other_side_crack_faces_nodes,
2085 other_side_crack_faces_nodes,
false);
2088 Range other_side_crack_faces_edges;
2091 other_side_crack_faces_edges,
2092 moab::Interface::UNION);
2093 other_side_crack_faces_edges =
2094 subtract(other_side_crack_faces_edges,
crackFront);
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);
2111 Range contact_both_sides_slave_nodes;
2114 Range contact_prisms;
2116 &*nit, 1, 3,
false, contact_prisms, moab::Interface::UNION);
2120 int side_number, other_side_number, sense, offset, number_nodes = 0;
2124 if (side_number < 3)
2125 contact_both_sides_slave_nodes.insert(conn[side_number + 3]);
2127 contact_both_sides_slave_nodes.insert(conn[side_number - 3]);
2129 contact_both_sides_slave_nodes =
2130 intersect(contact_both_sides_slave_nodes, prisms_level_nodes);
2133 CHKERR write_file(
"contact_both_sides_slave_nodes.vtk",
2134 contact_both_sides_slave_nodes);
2137 CHKERR add_both_sides_field(
"LAMBDA_BOTH_SIDES_CONTACT",
2139 contact_both_sides_slave_nodes,
false);
2146 const bool build_fields,
2156 Range crack_front_nodes;
2157 CHKERR moab.get_connectivity(crack_front_edges, crack_front_nodes,
true);
2160 1, MB_TAG_SPARSE,
MF_ZERO, verb);
2163 Range ents_to_remove;
2166 ents_to_remove = subtract(ents_to_remove, crack_front_nodes);
2168 ents_to_remove, verb);
2170 "LAMBDA_CRACKFRONT_AREA", verb);
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);
2186 Range ents_to_remove;
2189 ents_to_remove = subtract(ents_to_remove, field_ents);
2191 ents_to_remove, verb);
2193 field_ents, MBVERTEX,
"LAMBDA_CRACKFRONT_AREA_TANGENT", verb);
2207 const BitRefLevel mask2,
const bool add_forces,
const bool proc_only,
2213 bit1, mask1, MBTET, tets_level1);
2215 tets_level1 = intersect(
bitProcEnts, tets_level1);
2217 Range tets_not_level2;
2219 bit2, mask2, MBTET, tets_not_level2);
2220 tets_not_level2 = subtract(tets_level1, tets_not_level2);
2225 "SPATIAL_POSITION");
2227 "SPATIAL_POSITION");
2229 "SPATIAL_POSITION");
2231 "MESH_NODE_POSITIONS");
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);
2250 "SPATIAL_POSITION");
2252 "SPATIAL_POSITION");
2254 "SPATIAL_POSITION");
2256 "MESH_NODE_POSITIONS");
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);
2278 Range proc_tris, skin_tris, tets_level;
2280 moab::Interface::UNION);
2283 bit1, mask1, MBTET, tets_level);
2285 CHKERR skin.find_skin(0, tets_level,
false, skin_tris);
2286 skin_tris = intersect(skin_tris, proc_tris);
2291 "SPATIAL_POSITION");
2293 "SPATIAL_POSITION");
2295 "SPATIAL_POSITION");
2297 "MESH_NODE_POSITIONS");
2300 "SKIN",
"EIGEN_SPATIAL_POSITIONS");
2306 "MESH_NODE_POSITIONS");
2309 auto cn_value_ptr = boost::make_shared<double>(
cnValue);
2326 "MORTAR_CONTACT",
"SPATIAL_POSITION",
"LAMBDA_CONTACT",
2328 "EIGEN_SPATIAL_POSITIONS");
2332 Range all_slave_tris;
2337 "CONTACT_POST_PROC",
"SPATIAL_POSITION",
"LAMBDA_CONTACT",
2338 "MESH_NODE_POSITIONS", all_slave_tris);
2349 Range arc_fe_meshsets;
2352 arc_fe_meshsets.insert(fit->get()->getEnt());
2353 if (
auto ptr = fit->get()->getPartProcPtr())
2364 "LAMBDA_ARC_LENGTH");
2366 "LAMBDA_ARC_LENGTH");
2369 "LAMBDA_ARC_LENGTH");
2381 "ARC_LENGTH",
false);
2386 if (
auto ptr = fit->get()->getPartProcPtr())
2394 const bool proc_only) {
2400 bit, mask, MBTET, tets_level);
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);
2413 rval = moab.get_connectivity(tets_level, dofs_nodes,
true);
2415 Range lower_dim_ents = unite(dofs_ents, dofs_nodes);
2417 #ifdef __ANALITICAL_TRACTION__
2421 "SPATIAL_POSITION");
2423 "SPATIAL_POSITION");
2425 "SPATIAL_POSITION");
2427 "MESH_NODE_POSITIONS");
2430 if (meshset_manager_ptr->
checkMeshset(
"ANALITICAL_TRACTION")) {
2433 &cubit_meshset_ptr);
2437 tris = intersect(tris, lower_dim_ents);
2439 "ANALITICAL_TRACTION");
2443 #endif //__ANALITICAL_TRACTION__
2447 "SPATIAL_POSITION");
2449 "SPATIAL_POSITION");
2451 "SPATIAL_POSITION");
2453 "MESH_NODE_POSITIONS");
2464 moab::Interface::UNION);
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);
2485 const string block_set_force_name(
"FORCE");
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);
2494 for (
unsigned int ii = 0; ii < mydata.size(); ii++) {
2495 force[ii] = mydata[ii];
2497 if (force.size() != 3) {
2499 "Force vector not given");
2513 "SPATIAL_POSITION");
2515 "SPATIAL_POSITION");
2517 "SPATIAL_POSITION");
2519 "MESH_NODE_POSITIONS");
2525 tris = intersect(lower_dim_ents, tris);
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) {
2536 tris = intersect(lower_dim_ents, tris);
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);
2551 for (
unsigned int ii = 0; ii < mydata.size(); ii++) {
2552 pressure[ii] = mydata[ii];
2554 if (pressure.empty()) {
2556 "Pressure not given");
2569 #ifdef __ANALITICAL_DISPLACEMENT__
2573 if (meshset_manager_ptr->
checkMeshset(
"ANALITICAL_DISP")) {
2576 &cubit_meshset_ptr);
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);
2596 #endif //__ANALITICAL_DISPLACEMENT__
2603 const bool proc_only) {
2610 bit, mask, MBTET, tets_level);
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);
2622 rval = moab.get_connectivity(tets_level, dofs_nodes,
true);
2624 Range lower_dim_ents = unite(dofs_ents, dofs_nodes);
2628 "SPATIAL_POSITION");
2630 "SPATIAL_POSITION");
2632 "SPATIAL_POSITION");
2634 "MESH_NODE_POSITIONS");
2636 "MESH_NODE_POSITIONS");
2638 "MESH_NODE_POSITIONS");
2646 tris = intersect(lower_dim_ents, tris);
2647 ents_to_add.merge(tris);
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);
2671 bit, mask, MBTET, tets_level);
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);
2683 rval = moab.get_connectivity(tets_level, dofs_nodes,
true);
2685 Range lower_dim_ents = unite(dofs_ents, dofs_nodes);
2689 "SPATIAL_POSITION");
2691 "SPATIAL_POSITION");
2693 "SPATIAL_POSITION");
2695 "MESH_NODE_POSITIONS");
2697 "MESH_NODE_POSITIONS");
2699 "MESH_NODE_POSITIONS");
2707 tris = intersect(lower_dim_ents, tris);
2708 ents_to_add.merge(tris);
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) {
2718 tris = intersect(lower_dim_ents, tris);
2719 ents_to_add.merge(tris);
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);
2738 const bool proc_only) {
2745 bit, mask, MBTET, tets_level);
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);
2757 rval = moab.get_connectivity(tets_level, dofs_nodes,
true);
2759 Range lower_dim_ents = unite(dofs_ents, dofs_nodes);
2763 "SPATIAL_POSITION");
2765 "SPATIAL_POSITION");
2767 "SPATIAL_POSITION");
2769 "MESH_NODE_POSITIONS");
2771 "MESH_NODE_POSITIONS");
2773 "MESH_NODE_POSITIONS");
2777 if (
bit->getName().compare(0, 9,
"SPRING_BC") == 0) {
2781 tris = intersect(lower_dim_ents, tris);
2782 ents_to_add.merge(tris);
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);
2807 bit, mask, MBPRISM, prisms_level);
2810 prisms_level = intersect(
bitProcEnts, prisms_level);
2816 "SIMPLE_CONTACT_ALE",
"SPATIAL_POSITION",
"LAMBDA_CONTACT",
2818 "EIGEN_SPATIAL_POSITIONS");
2822 Range face_on_prism;
2823 CHKERR moab.get_adjacencies(contact_prisms, 2,
false, face_on_prism,
2824 moab::Interface::UNION);
2826 Range tris_on_prism = face_on_prism.subset_by_type(MBTRI);
2831 Range adj_vols_to_prisms;
2832 CHKERR moab.get_adjacencies(check_tris, 3,
false, adj_vols_to_prisms,
2833 moab::Interface::UNION);
2837 bit, mask, MBTET, tet_level);
2839 Range contact_tets_from_vols = adj_vols_to_prisms.subset_by_type(MBTET);
2841 Range contact_tets = intersect(tet_level, contact_tets_from_vols);
2845 "SPATIAL_POSITION");
2847 "SPATIAL_POSITION");
2849 "MESH_NODE_POSITIONS");
2851 "MESH_NODE_POSITIONS");
2853 "SPATIAL_POSITION");
2855 "MESH_NODE_POSITIONS");
2857 Range 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);
2874 const bool proc_only,
2880 bit, mask, MBTET, tets_level);
2888 "SPATIAL_POSITION");
2890 "SPATIAL_POSITION");
2892 "MESH_NODE_POSITIONS");
2894 "MESH_NODE_POSITIONS");
2896 "SPATIAL_POSITION");
2898 "MESH_NODE_POSITIONS");
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);
2919 const bool proc_only,
2926 "MESH_NODE_POSITIONS");
2928 "MESH_NODE_POSITIONS");
2930 "MESH_NODE_POSITIONS");
2931 Range crack_front_adj_tris;
2933 crackFrontNodes, 2,
false, crack_front_adj_tris, moab::Interface::UNION);
2934 crack_front_adj_tris = intersect(crack_front_adj_tris,
crackFaces);
2937 bit, mask, MBTRI, tris_level);
2938 crack_front_adj_tris = intersect(crack_front_adj_tris, tris_level);
2940 Range ents_to_remove;
2943 ents_to_remove = subtract(ents_to_remove, crack_front_adj_tris);
2945 ents_to_remove, verb);
2948 "GRIFFITH_FORCE_ELEMENT");
2950 &crack_front_adj_tris, verb);
2952 #ifdef __ANALITICAL_TRACTION__
2956 bit, mask, MBTET, tets_level);
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);
2968 CHKERR moab.get_connectivity(tets_level, dofs_nodes,
true);
2969 Range lower_dim_ents = unite(dofs_ents, dofs_nodes);
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");
2983 if (meshset_manager_ptr->
checkMeshset(
"ANALITICAL_TRACTION")) {
2986 &cubit_meshset_ptr);
2990 tris = intersect(tris, lower_dim_ents);
2992 tris, MBTRI,
"ANALITICAL_METERIAL_TRACTION");
2997 #endif //__ANALITICAL_TRACTION__
3003 const BitRefLevel mask,
const bool proc_only,
const bool verb) {
3006 auto get_prisms_level = [&](
auto bit) {
3009 bit, mask, MBPRISM, prisms_level);
3011 prisms_level = intersect(
bitProcEnts, prisms_level);
3013 return prisms_level;
3016 auto get_crack_prisms = [&](
auto prisms_level) {
3019 return crack_prisms;
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) {
3029 if (!ents.empty()) {
3052 CHKERR add_both_sides_fe(
"CLOSE_CRACK",
"LAMBDA_CLOSE_CRACK",
3054 get_crack_prisms(get_prisms_level(bit_spatial)));
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)));
3060 Range contact_prisms =
3062 CHKERR add_both_sides_fe(
"BOTH_SIDES_OF_CONTACT",
3063 "LAMBDA_BOTH_SIDES_CONTACT",
"MESH_NODE_POSITIONS",
3071 const bool proc_only,
3076 bit, mask, MBTET, tets_level);
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");
3097 Range ents_to_remove;
3100 ents_to_remove = subtract(ents_to_remove, tets_level);
3112 const std::vector<int> &ids,
const bool proc_only,
const int verb,
3130 moab::Interface::UNION);
3131 surface_ents = intersect(surface_ents, proc_tris);
3136 "MESH_NODE_POSITIONS");
3138 "MESH_NODE_POSITIONS");
3140 "MESH_NODE_POSITIONS");
3142 for (std::vector<int>::const_iterator it = ids.begin(); it != ids.end();
3145 "LAMBDA_SURFACE" + boost::lexical_cast<std::string>(*it);
3150 CHKERR moab.get_adjacencies(&ent, 1, 2,
false, adj_tris);
3151 adj_tris = intersect(adj_tris,
bodySkin);
3152 surface_ents.merge(adj_tris);
3164 Range ents_to_remove;
3167 ents_to_remove = subtract(ents_to_remove, surface_ents);
3181 const int verb,
const bool debug) {
3186 auto get_edges_block_ents = [
this, meshset_manager_ptr]() {
3187 EntityHandle edge_block_set = getInterface<CPMeshCut>()->getEdgesBlockSet();
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);
3193 return meshset_ents;
3196 auto get_skin = [
this, bit_ref_manager_ptr]() {
3198 CHKERR bit_ref_manager_ptr->getEntitiesByTypeAndRefLevel(
3202 CHKERR skin.find_skin(0, tets,
false, skin_faces);
3206 auto intersect_with_skin_edges = [
this, get_edges_block_ents, get_skin]() {
3209 moab::Interface::UNION);
3210 return intersect(get_edges_block_ents(), adj_edges);
3212 Range edges_ents = intersect_with_skin_edges();
3214 Range contact_faces;
3219 boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>(
3229 Range faces = get_skin();
3234 auto intersect_with_bit_level_edges = [bit_ref_manager_ptr, &
bit,
3237 CHKERR bit_ref_manager_ptr->getEntitiesByTypeAndRefLevel(
3239 return intersect(edges_ents, level_edges);
3241 edges_ents = intersect_with_bit_level_edges();
3243 auto intersect_with_proc_edges = [
this, &edges_ents, proc_only]() {
3247 bitProcEnts, 1,
false, proc_edges, moab::Interface::UNION);
3248 return intersect(edges_ents, proc_edges);
3252 edges_ents = intersect_with_proc_edges();
3259 auto add_field = [
this, &fe_name](std::string &&
field_name) {
3266 CHKERR add_field(
"MESH_NODE_POSITIONS");
3267 CHKERR add_field(
"LAMBDA_EDGE");
3272 "edges_normals_on_proc_" +
3278 Range ents_to_remove;
3280 ents_to_remove = subtract(ents_to_remove, edges_ents);
3291 const bool proc_only,
const int verb,
const bool debug) {
3299 bit, mask, MBTRI, level_tris);
3300 Range crack_surface;
3305 moab::Interface::UNION);
3306 crack_surface = intersect(crack_surface, proc_tris);
3308 Range crack_front_adj_tris;
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);
3315 "LAMBDA_SURFACE" + boost::lexical_cast<std::string>(
3316 getInterface<CPMeshCut>()->getCrackSurfaceId());
3323 "MESH_NODE_POSITIONS");
3327 "MESH_NODE_POSITIONS");
3331 "MESH_NODE_POSITIONS");
3334 crack_surface,
QUIET);
3340 ents = subtract(ents, crack_surface);
3341 if (!ents.empty()) {
3351 std::ostringstream file;
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);
3363 "LAMBDA_CRACKFRONT_AREA");
3365 "MESH_NODE_POSITIONS");
3367 "LAMBDA_CRACKFRONT_AREA");
3369 "MESH_NODE_POSITIONS");
3371 "MESH_NODE_POSITIONS");
3373 "LAMBDA_CRACKFRONT_AREA");
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");
3392 crack_front_adj_tris,
QUIET);
3394 auto remeve_ents_from_fe =
3395 [
this, &crack_front_adj_tris](
const std::string fe_name) {
3399 ents = subtract(ents, crack_front_adj_tris);
3400 if (!ents.empty()) {
3405 CHKERR remeve_ents_from_fe(
"CRACKFRONT_AREA_ELEM");
3406 CHKERR remeve_ents_from_fe(
"CRACKFRONT_AREA_TANGENT_ELEM");
3409 "CRACKFRONT_AREA_ELEM");
3411 crack_front_adj_tris, MBTRI,
"CRACKFRONT_AREA_TANGENT_ELEM");
3414 std::ostringstream file;
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);
3424 &crack_front_adj_tris, verb);
3426 &crack_front_adj_tris, verb);
3435 const std::vector<std::string> fe_list) {
3471 for (std::vector<std::string>::const_iterator it = fe_list.begin();
3472 it != fe_list.end(); it++)
3478 #ifdef __ANALITICAL_TRACTION__
3500 PetscSection section;
3502 CHKERR DMSetDefaultSection(dm, section);
3503 CHKERR DMSetDefaultGlobalSection(dm, section);
3505 CHKERR PetscSectionDestroy(§ion);
3511 const std::string prb_name,
3529 #ifdef __ANALITICAL_TRACTION__
3540 PetscSection section;
3542 CHKERR DMSetDefaultSection(dm, section);
3543 CHKERR DMSetDefaultGlobalSection(dm, section);
3545 CHKERR PetscSectionDestroy(§ion);
3552 arcCtx = boost::make_shared<ArcLengthCtx>(
mField, problem_ptr->getName(),
3553 "LAMBDA_ARC_LENGTH");
3555 auto arc_snes_ctx = boost::make_shared<ArcLengthSnesCtx>(
3581 #ifdef __ANALITICAL_TRACTION__
3590 auto remove_higher_dofs_than_approx_dofs_from_eigen_elastic = [&]() {
3595 CHKERR prb_mng->removeDofsOnEntities(
"EIGEN_ELASTIC",
"SPATIAL_POSITION",
3599 CHKERR remove_higher_dofs_than_approx_dofs_from_eigen_elastic();
3603 PetscSection section;
3605 CHKERR DMSetDefaultSection(dm, section);
3606 CHKERR DMSetDefaultGlobalSection(dm, section);
3608 CHKERR PetscSectionDestroy(§ion);
3616 "LAMBDA_ARC_LENGTH");
3618 auto arc_snes_ctx = boost::make_shared<ArcLengthSnesCtx>(
3626 const std::string prb_name,
3635 #ifdef __ANALITICAL_DISPLACEMENT__
3636 if (meshset_manager_ptr->
checkMeshset(
"ANALITICAL_DISP")) {
3639 &cubit_meshset_ptr);
3643 if (!tris.empty()) {
3655 #endif //__ANALITICAL_DISPLACEMENT__
3661 const BitRefLevel mask,
const std::vector<std::string> fe_list,
3669 for (
auto &fe : fe_list)
3675 #ifdef __ANALITICAL_TRACTION__
3694 PetscSection section;
3696 CHKERR DMSetDefaultSection(dm, section);
3697 CHKERR DMSetDefaultGlobalSection(dm, section);
3699 CHKERR PetscSectionDestroy(§ion);
3702 std::vector<std::string> fe_lists;
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();
3713 prb_name, fe_name, &meshset);
3715 (prb_name +
"_" + fe_name +
".h5m").c_str(),
"MOAB",
3716 "PARALLEL=WRITE_PART", &meshset, 1);
3726 const std::string prb_name,
const int verb) {
3738 #ifdef __ANALITICAL_TRACTION__
3749 const std::string prb_name,
const std::vector<int> surface_ids,
3750 std::vector<std::string> fe_list,
const int verb) {
3756 for (
auto id : surface_ids) {
3758 dm, (
"LAMBDA_SURFACE" + boost::lexical_cast<std::string>(
id)).c_str());
3761 for (
auto fe : fe_list) {
3764 auto get_edges_block_set = [
this]() {
3765 return getInterface<CPMeshCut>()->getEdgesBlockSet();
3767 if (get_edges_block_set()) {
3778 const std::string prb_name,
const bool verb) {
3796 const int verb,
const bool debug) {
3804 feRhs = boost::make_shared<CrackFrontElement>(
3809 feLhs = boost::make_shared<CrackFrontElement>(
3815 feRhs->meshPositionsFieldName =
"NONE";
3816 feLhs->meshPositionsFieldName =
"NONE";
3817 feRhs->addToRule = 0;
3818 feLhs->addToRule = 0;
3823 boost::shared_ptr<map<int, BlockData>> block_sets_ptr(
3831 CHKERR it->getAttributeDataStructure(mydata);
3832 int id = it->getMeshsetId();
3835 meshset, MBTET,
elasticFe->setOfBlocks[
id].tEts,
true);
3838 elasticFe->setOfBlocks[id].PoissonRatio = mydata.
data.Poisson;
3841 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"\nMaterial block %d \n",
id);
3842 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"\tYoung's modulus %6.4e \n",
3844 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"\tPoisson's ratio %6.4e \n\n",
3845 mydata.
data.Poisson);
3850 elasticFe->setOfBlocks[id].materialDoublePtr =
3851 boost::make_shared<Hooke<double>>();
3852 elasticFe->setOfBlocks[id].materialAdoublePtr =
3853 boost::make_shared<Hooke<adouble>>();
3856 elasticFe->setOfBlocks[id].materialDoublePtr = boost::make_shared<
3859 elasticFe->setOfBlocks[id].materialAdoublePtr = boost::make_shared<
3865 elasticFe->setOfBlocks[id].materialDoublePtr =
3866 boost::make_shared<NeoHookean<double>>();
3867 elasticFe->setOfBlocks[id].materialAdoublePtr =
3868 boost::make_shared<NeoHookean<adouble>>();
3872 elasticFe->setOfBlocks[id].materialDoublePtr =
3873 boost::make_shared<Hooke<double>>();
3874 elasticFe->setOfBlocks[id].materialAdoublePtr =
3875 boost::make_shared<Hooke<adouble>>();
3879 "Material type not yet implemented");
3886 elasticFe->commonData.spatialPositions =
"SPATIAL_POSITION";
3887 elasticFe->commonData.meshPositions =
"MESH_NODE_POSITIONS";
3889 auto data_hooke_element_at_pts =
3890 boost::make_shared<HookeElement::DataAtIntegrationPts>();
3893 auto mat_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
3895 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
3896 boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
3898 mat_grad_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
3899 feRhs->getOpPtrVector().push_back(
3901 mat_grad_pos_at_pts_ptr));
3910 "SPATIAL_POSITION",
elasticFe->commonData,
feRhs->singularElement,
3913 feRhs->getOpPtrVector().push_back(
3915 "MESH_NODE_POSITIONS",
elasticFe->commonData));
3920 map<int, NonlinearElasticElement::BlockData>::iterator sit =
3922 for (; sit !=
elasticFe->setOfBlocks.end(); sit++) {
3924 feRhs->getOpPtrVector().push_back(
3926 "SPATIAL_POSITION", sit->second,
elasticFe->commonData,
3929 feRhs->getOpPtrVector().push_back(
3931 "SPATIAL_POSITION", sit->second,
elasticFe->commonData));
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;
3952 th_rho) != MB_SUCCESS) {
3954 "Density tag name %s cannot be found. Please "
3955 "provide the correct name.",
3969 feRhs->getOpPtrVector().push_back(
3971 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
3974 if (
mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
3979 feRhs->getOpPtrVector().push_back(
3981 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
feRhs,
3983 feRhs->getOpPtrVector().push_back(
3985 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
feRhs,
3989 feRhs->getOpPtrVector().push_back(
3990 new HookeElement::OpCalculateStiffnessScaledByDensityField(
3991 "SPATIAL_POSITION",
"SPATIAL_POSITION", block_sets_ptr,
3994 feRhs->getOpPtrVector().push_back(
3995 new HookeElement::OpCalculateStrainAle(
"SPATIAL_POSITION",
3997 data_hooke_element_at_pts));
3998 feRhs->getOpPtrVector().push_back(
3999 new HookeElement::OpCalculateStress<1>(
"SPATIAL_POSITION",
4001 data_hooke_element_at_pts));
4003 feRhs->getOpPtrVector().push_back(
4004 new HookeElement::OpCalculateHomogeneousStiffness<0>(
4005 "SPATIAL_POSITION",
"SPATIAL_POSITION", block_sets_ptr,
4006 data_hooke_element_at_pts));
4010 feRhs->getOpPtrVector().push_back(
4012 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
4014 feRhs->getOpPtrVector().push_back(
4015 new HookeElement::OpCalculateStrainAle(
"SPATIAL_POSITION",
4017 data_hooke_element_at_pts));
4018 feRhs->getOpPtrVector().push_back(
4019 new HookeElement::OpCalculateStress<0>(
"SPATIAL_POSITION",
4021 data_hooke_element_at_pts));
4026 feRhs->getOpPtrVector().push_back(
new HookeElement::OpAleRhs_dx(
4027 "SPATIAL_POSITION",
"SPATIAL_POSITION", data_hooke_element_at_pts));
4042 if (
mwlsApprox->mwlsMoab.tag_get_handle(mwls_stress_tag_name.c_str(),
4043 th) != MB_SUCCESS) {
4045 "Internal stress tag name %s cannot be found. Please "
4046 "provide the correct name.",
4047 mwls_stress_tag_name.substr(4).c_str());
4052 if (
mwlsApprox->mwlsMoab.tag_get_handle(mwls_stress_tag_name.c_str(),
4053 th) != MB_SUCCESS) {
4055 "Internal stress tag name %s cannot be found. Please "
4056 "provide the correct name.",
4057 mwls_stress_tag_name.substr(4).c_str());
4064 if (
mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
4065 feRhs->getOpPtrVector().push_back(
4068 mwls_stress_tag_name,
false));
4070 feRhs->getOpPtrVector().push_back(
4072 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
feRhs,
4074 feRhs->getOpPtrVector().push_back(
4077 mwls_stress_tag_name,
false));
4084 feRhs->getOpPtrVector().push_back(
4085 new HookeElement::OpAnalyticalInternalAleStrain_dx<0>(
4086 "SPATIAL_POSITION", data_hooke_element_at_pts,
4088 boost::shared_ptr<MatrixDouble>(
4099 "SPATIAL_POSITION",
elasticFe->commonData,
feLhs->singularElement,
4102 feLhs->getOpPtrVector().push_back(
4104 "MESH_NODE_POSITIONS",
elasticFe->commonData));
4109 map<int, NonlinearElasticElement::BlockData>::iterator sit =
4111 for (; sit !=
elasticFe->setOfBlocks.end(); sit++) {
4113 feLhs->getOpPtrVector().push_back(
4115 "SPATIAL_POSITION", sit->second,
elasticFe->commonData,
4118 feLhs->getOpPtrVector().push_back(
4120 "SPATIAL_POSITION",
"SPATIAL_POSITION", sit->second,
4126 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
4127 boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
4129 mat_grad_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
4130 feLhs->getOpPtrVector().push_back(
4132 mat_grad_pos_at_pts_ptr));
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;
4143 "mwlsApprox not allocated");
4154 feLhs->getOpPtrVector().push_back(
4156 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
4160 if (
mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
4165 feLhs->getOpPtrVector().push_back(
4167 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
feLhs,
4169 feLhs->getOpPtrVector().push_back(
4171 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
feLhs,
4175 feLhs->getOpPtrVector().push_back(
4176 new HookeElement::OpCalculateStiffnessScaledByDensityField(
4177 "SPATIAL_POSITION",
"SPATIAL_POSITION", block_sets_ptr,
4180 feLhs->getOpPtrVector().push_back(
4181 new HookeElement::OpCalculateStrainAle(
"SPATIAL_POSITION",
4183 data_hooke_element_at_pts));
4184 feLhs->getOpPtrVector().push_back(
4185 new HookeElement::OpCalculateStress<1>(
"SPATIAL_POSITION",
4187 data_hooke_element_at_pts));
4191 feLhs->getOpPtrVector().push_back(
new HookeElement::OpAleLhs_dx_dx<1>(
4192 "SPATIAL_POSITION",
"SPATIAL_POSITION", data_hooke_element_at_pts));
4195 feLhs->getOpPtrVector().push_back(
4196 new HookeElement::OpCalculateHomogeneousStiffness<0>(
4197 "SPATIAL_POSITION",
"SPATIAL_POSITION", block_sets_ptr,
4198 data_hooke_element_at_pts));
4202 feLhs->getOpPtrVector().push_back(
4204 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
4206 feLhs->getOpPtrVector().push_back(
4207 new HookeElement::OpCalculateStrainAle(
"SPATIAL_POSITION",
4209 data_hooke_element_at_pts));
4210 feLhs->getOpPtrVector().push_back(
4211 new HookeElement::OpCalculateStress<0>(
"SPATIAL_POSITION",
4213 data_hooke_element_at_pts));
4218 feLhs->getOpPtrVector().push_back(
new HookeElement::OpAleLhs_dx_dx<0>(
4219 "SPATIAL_POSITION",
"SPATIAL_POSITION", data_hooke_element_at_pts));
4227 boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
4230 string fe_name_str =
"FORCE_FE";
4233 auto &fe = nf.getLoopFe();
4240 .addForce(
"SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(),
true);
4243 const string block_set_force_name(
"FORCE");
4246 if (it->getName().compare(0, block_set_force_name.length(),
4247 block_set_force_name) == 0) {
4249 .addForce(
"SPATIAL_POSITION", PETSC_NULL, (it->getMeshsetId()),
4257 boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
4259 boost::make_shared<NeumannForcesSurface::DataAtIntegrationPts>();
4262 string fe_name_str =
"FORCE_FE_ALE";
4266 auto &fe_lhs = nf.getLoopFeLhs();
4267 auto &fe_mat_rhs = nf.getLoopFeMatRhs();
4268 auto &fe_mat_lhs = nf.getLoopFeMatLhs();
4277 .addForceAle(
"SPATIAL_POSITION",
"MESH_NODE_POSITIONS",
4279 PETSC_NULL,
bit->getMeshsetId(),
true,
false,
4283 const string block_set_force_name(
"FORCE");
4286 if (it->getName().compare(0, block_set_force_name.length(),
4287 block_set_force_name) == 0) {
4289 .addForceAle(
"SPATIAL_POSITION",
"MESH_NODE_POSITIONS",
4291 PETSC_NULL, it->getMeshsetId(),
true,
true,
4299 boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
4301 string fe_name_str =
"PRESSURE_FE";
4305 auto &fe = nf.getLoopFe();
4313 .addPressure(
"SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(),
4317 const string block_set_pressure_name(
"PRESSURE");
4319 if (it->getName().compare(0, block_set_pressure_name.length(),
4320 block_set_pressure_name) == 0) {
4322 .addPressure(
"SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(),
4327 const string block_set_linear_pressure_name(
"LINEAR_PRESSURE");
4329 if (it->getName().compare(0, block_set_linear_pressure_name.length(),
4330 block_set_linear_pressure_name) == 0) {
4332 .addLinearPressure(
"SPATIAL_POSITION", PETSC_NULL,
4333 it->getMeshsetId(),
true);
4340 boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
4342 boost::make_shared<NeumannForcesSurface::DataAtIntegrationPts>();
4344 string fe_name_str =
"PRESSURE_ALE";
4348 auto &fe_lhs = nf.getLoopFeLhs();
4349 auto &fe_mat_rhs = nf.getLoopFeMatRhs();
4350 auto &fe_mat_lhs = nf.getLoopFeMatLhs();
4360 .addPressureAle(
"SPATIAL_POSITION",
"MESH_NODE_POSITIONS",
4362 PETSC_NULL, it->getMeshsetId(),
true);
4368 boost::make_shared<FaceElementForcesAndSourcesCore>(
mField);
4371 boost::make_shared<FaceElementForcesAndSourcesCore>(
mField);
4373 boost::make_shared<FaceElementForcesAndSourcesCore>(
mField);
4375 boost::make_shared<MetaSpringBC::DataAtIntegrationPtsSprings>(
mField);
4386 feSpringLhsPtr = boost::shared_ptr<FaceElementForcesAndSourcesCore>(
4388 feSpringRhsPtr = boost::shared_ptr<FaceElementForcesAndSourcesCore>(
4393 "MESH_NODE_POSITIONS");
4395 bool use_eigen_pos =
4398 bool use_eigen_pos_simple_contact =
4403 boost::make_shared<SimpleContactProblem::SimpleContactElement>(
mField);
4405 boost::make_shared<SimpleContactProblem::SimpleContactElement>(
mField);
4407 boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(
mField);
4416 "LAMBDA_CONTACT",
false, use_eigen_pos_simple_contact,
4417 "EIGEN_SPATIAL_POSITIONS");
4420 "LAMBDA_CONTACT",
false, use_eigen_pos_simple_contact,
4421 "EIGEN_SPATIAL_POSITIONS");
4425 "LAMBDA_CONTACT",
false, use_eigen_pos_simple_contact,
4426 "EIGEN_SPATIAL_POSITIONS");
4430 "LAMBDA_CONTACT",
false, use_eigen_pos_simple_contact,
4431 "EIGEN_SPATIAL_POSITIONS");
4438 mField,
"LAMBDA_CLOSE_CRACK",
"SPATIAL_POSITION"));
4447 boost::make_shared<SimpleContactProblem::SimpleContactElement>(
mField);
4449 boost::make_shared<SimpleContactProblem::SimpleContactElement>(
mField);
4451 boost::make_shared<SimpleContactProblem::SimpleContactElement>(
mField);
4453 boost::make_shared<SimpleContactProblem::CommonDataSimpleContact>(
4460 "SPATIAL_POSITION",
"MESH_NODE_POSITIONS",
"LAMBDA_CONTACT",
4465 "SPATIAL_POSITION",
"MESH_NODE_POSITIONS",
"LAMBDA_CONTACT",
4470 "MESH_NODE_POSITIONS",
"LAMBDA_CONTACT", use_eigen_pos_simple_contact,
4471 "EIGEN_SPATIAL_POSITIONS");
4475 boost::make_shared<MortarContactProblem::MortarContactElement>(
4477 "MESH_NODE_POSITIONS");
4479 boost::make_shared<MortarContactProblem::MortarContactElement>(
4481 "MESH_NODE_POSITIONS");
4483 boost::make_shared<MortarContactProblem::CommonDataMortarContact>(
mField);
4492 "LAMBDA_CONTACT",
false, use_eigen_pos,
"EIGEN_SPATIAL_POSITIONS");
4496 "LAMBDA_CONTACT",
false, use_eigen_pos,
"EIGEN_SPATIAL_POSITIONS");
4500 "LAMBDA_CONTACT",
false, use_eigen_pos,
"EIGEN_SPATIAL_POSITIONS");
4504 "LAMBDA_CONTACT",
false, use_eigen_pos,
"EIGEN_SPATIAL_POSITIONS");
4510 boost::make_shared<SimpleContactProblem::SimpleContactElement>(
mField);
4513 boost::make_shared<MortarContactProblem::MortarContactElement>(
4515 "MESH_NODE_POSITIONS");
4521 use_eigen_pos_simple_contact,
"EIGEN_SPATIAL_POSITIONS");
4526 use_eigen_pos,
"EIGEN_SPATIAL_POSITIONS");
4530 nodalForces = boost::make_shared<boost::ptr_map<string, NodalForce>>();
4532 string fe_name_str =
"FORCE_FE";
4537 .addForce(
"SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(),
false);
4541 edgeForces = boost::make_shared<boost::ptr_map<string, EdgeForce>>();
4543 string fe_name_str =
"FORCE_FE";
4548 .addForce(
"SPATIAL_POSITION", PETSC_NULL, it->getMeshsetId(),
false);
4555 PETSC_NULL, PETSC_NULL));
4558 #ifdef __ANALITICAL_DISPLACEMENT__
4561 mField,
"SPATIAL_POSITION", PETSC_NULL, PETSC_NULL, PETSC_NULL));
4562 #endif //__ANALITICAL_DISPLACEMENT__
4565 #ifdef __ANALITICAL_TRACTION__
4568 if (meshset_manager_ptr->
checkMeshset(
"ANALITICAL_TRACTION")) {
4571 boost::make_shared<boost::ptr_map<string, NeumannForcesSurface>>();
4573 string fe_name_str =
"ANALITICAL_TRACTION";
4579 &cubit_meshset_ptr);
4586 .fe.getOpPtrVector()
4590 "SPATIAL_POSITION", faces,
4596 #endif //__ANALITICAL_TRACTION__
4611 DM dm, Mat
m,
Vec q,
Vec f, boost::shared_ptr<FEMethod> arc_method,
4612 boost::shared_ptr<ArcLengthCtx> arc_ctx,
const int verb,
const bool debug) {
4613 boost::shared_ptr<FEMethod>
null;
4629 VecSetOption(arc_ctx->F_lambda, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
4634 auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
4635 auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
4636 for (; oit != hi_oit; oit++) {
4637 if (boost::typeindex::type_id_runtime(*oit) ==
4638 boost::typeindex::type_id<NeumannForcesSurface::OpNeumannForce>()) {
4646 auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
4647 auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
4648 for (; oit != hi_oit; oit++) {
4649 if (boost::typeindex::type_id_runtime(*oit) ==
4650 boost::typeindex::type_id<
4655 if (boost::typeindex::type_id_runtime(*oit) ==
4656 boost::typeindex::type_id<
4664 auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
4665 auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
4666 for (; oit != hi_oit; oit++) {
4667 if (boost::typeindex::type_id_runtime(*oit) ==
4668 boost::typeindex::type_id<EdgeForce::OpEdgeForce>()) {
4674 auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
4675 auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
4676 for (; oit != hi_oit; oit++) {
4677 if (boost::typeindex::type_id_runtime(*oit) ==
4678 boost::typeindex::type_id<NodalForce::OpNodalForce>()) {
4687 #ifdef __ANALITICAL_DISPLACEMENT__
4692 #endif //__ANALITICAL_DISPLACEMENT__
4701 if (problem_ptr->
getName() ==
"EIGEN_ELASTIC") {
4707 #ifdef __ANALITICAL_DISPLACEMENT__
4712 #endif //__ANALITICAL_DISPLACEMENT__
4718 auto fe_set_option = boost::make_shared<FEMethod>();
4719 fe_set_option->preProcessHook = [fe_set_option]() {
4720 return VecSetOption(fe_set_option->snes_f, VEC_IGNORE_NEGATIVE_INDICES,
4727 #ifdef __ANALITICAL_DISPLACEMENT__
4732 #endif //__ANALITICAL_DISPLACEMENT__
4741 if (problem_ptr->
getName() ==
"EIGEN_ELASTIC") {
4749 dm, fit->first.c_str(),
4750 boost::shared_ptr<FEMethod>(
surfaceForces, &fit->second->getLoopFe()),
4756 dm, fit->first.c_str(),
4757 boost::shared_ptr<FEMethod>(
surfacePressure, &fit->second->getLoopFe()),
4762 dm, fit->first.c_str(),
4763 boost::shared_ptr<FEMethod>(
edgeForces, &fit->second->getLoopFe()),
4768 dm, fit->first.c_str(),
4769 boost::shared_ptr<FEMethod>(
nodalForces, &fit->second->getLoopFe()),
4772 #ifdef __ANALITICAL_TRACTION__
4777 dm, fit->first.c_str(),
4783 #endif //__ANALITICAL_TRACTION__
4784 #ifdef __ANALITICAL_DISPLACEMENT__
4789 #endif //__ANALITICAL_DISPLACEMENT__
4796 if (
m == PETSC_NULL || q == PETSC_NULL ||
f == PETSC_NULL) {
4798 "problem matrix or vectors are set to null");
4801 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
4808 CHKERR SNESDestroy(&snes);
4809 MatView(
m, PETSC_VIEWER_DRAW_WORLD);
4822 boost::shared_ptr<NonlinearElasticElement>
fE;
4832 if (
type != MBVERTEX)
4854 cerr << data.getDiffN() << endl;
4860 gg !=
static_cast<int>(
4861 fE->commonData.gradAtGaussPts[
"MESH_NODE_POSITIONS"].size());
4864 fE->commonData.gradAtGaussPts[
"MESH_NODE_POSITIONS"][gg];
4866 &mat(1, 0), &mat(1, 1), &mat(1, 2),
4867 &mat(2, 0), &mat(2, 1), &mat(2, 2));
4871 cerr <<
"[" << tH(0, 0) <<
"," << tH(0, 1) <<
"," << tH(0, 2) <<
";"
4872 << tH(1, 0) <<
"," << tH(1, 1) <<
"," << tH(1, 2) <<
";" << tH(2, 0)
4873 <<
"," << tH(2, 1) <<
"," << tH(2, 2) <<
"]" << endl;
4903 boost::shared_ptr<CrackFrontElement> &fe_rhs,
4904 boost::shared_ptr<CrackFrontElement> &fe_lhs) {
4912 "Elastic element instance not created");
4917 fe_rhs = boost::make_shared<CrackFrontElement>(
4921 fe_lhs = boost::make_shared<CrackFrontElement>(
4927 fe_rhs->meshPositionsFieldName =
"NONE";
4928 fe_lhs->meshPositionsFieldName =
"NONE";
4929 fe_rhs->addToRule = 0;
4930 fe_lhs->addToRule = 0;
4945 "Internal stress tag name %s cannot be found. Please "
4946 "provide the correct name.",
4957 auto mat_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
4959 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
4961 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
4962 boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr(
new MatrixDouble());
4964 "MESH_NODE_POSITIONS", mat_grad_pos_at_pts_ptr));
4966 "MESH_NODE_POSITIONS", mat_grad_pos_at_pts_ptr));
4968 boost::shared_ptr<MatrixDouble> space_grad_pos_at_pts_ptr(
new MatrixDouble());
4970 "SPATIAL_POSITION", space_grad_pos_at_pts_ptr, fe_rhs->singularElement,
4973 fe_rhs->singularElement, fe_rhs->detS, fe_rhs->invSJac));
4975 if (
mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
4977 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_rhs,
mwlsApprox,
4980 fe_rhs->getOpPtrVector().push_back(
4982 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_rhs,
mwlsApprox));
4983 fe_rhs->getOpPtrVector().push_back(
4985 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_rhs,
mwlsApprox,
4990 mat_grad_pos_at_pts_ptr,
mwlsApprox,
false));
4992 space_grad_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
mwlsApprox,
4996 "SPATIAL_POSITION", space_grad_pos_at_pts_ptr, fe_lhs->singularElement,
4999 fe_lhs->singularElement, fe_lhs->detS, fe_lhs->invSJac));
5001 if (
mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5003 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_lhs,
mwlsApprox,
5006 fe_lhs->getOpPtrVector().push_back(
5008 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_lhs,
mwlsApprox));
5009 fe_lhs->getOpPtrVector().push_back(
5011 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_lhs,
mwlsApprox,
5018 space_grad_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
mwlsApprox,
5021 space_grad_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
mwlsApprox,
5028 boost::shared_ptr<CrackFrontElement> &fe_rhs,
5029 boost::shared_ptr<CrackFrontElement> &fe_lhs) {
5037 "Elastic element instance not created");
5042 fe_rhs = boost::make_shared<CrackFrontElement>(
5046 fe_lhs = boost::make_shared<CrackFrontElement>(
5050 bool test_with_mwls =
false;
5054 fe_rhs->meshPositionsFieldName =
"NONE";
5055 fe_lhs->meshPositionsFieldName =
"NONE";
5056 fe_rhs->addToRule = 0;
5057 fe_lhs->addToRule = 0;
5059 boost::shared_ptr<HookeElement::DataAtIntegrationPts>
5060 data_hooke_element_at_pts(
new HookeElement::DataAtIntegrationPts());
5061 boost::shared_ptr<map<int, NonlinearElasticElement::BlockData>>
5077 "Density tag name %s cannot be found. Please "
5078 "provide the correct name.",
5089 auto mat_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
5091 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
5093 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
5095 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_MWLS_ptr(
5098 "MESH_NODE_POSITIONS", mat_pos_at_pts_MWLS_ptr));
5100 "MESH_NODE_POSITIONS", mat_pos_at_pts_MWLS_ptr));
5101 auto mat_grad_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
5103 "MESH_NODE_POSITIONS", mat_grad_pos_at_pts_ptr));
5105 "MESH_NODE_POSITIONS", mat_grad_pos_at_pts_ptr));
5107 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
5109 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
5111 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
5112 fe_rhs->singularElement, fe_rhs->invSJac));
5114 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
5115 fe_lhs->singularElement, fe_lhs->invSJac));
5116 auto space_grad_pos_at_pts_ptr = boost::make_shared<MatrixDouble>();
5118 "SPATIAL_POSITION", space_grad_pos_at_pts_ptr, fe_rhs->singularElement,
5121 "SPATIAL_POSITION", space_grad_pos_at_pts_ptr, fe_lhs->singularElement,
5125 fe_rhs->singularElement, fe_rhs->detS, fe_rhs->invSJac));
5131 if (test_with_mwls) {
5132 if (
mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5134 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_rhs,
mwlsApprox,
5137 fe_rhs->getOpPtrVector().push_back(
5139 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_rhs,
mwlsApprox));
5140 fe_rhs->getOpPtrVector().push_back(
5142 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_rhs,
mwlsApprox,
5148 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr,
mwlsApprox->rhoAtGaussPts,
5150 fe_rhs,
mwlsApprox->singularInitialDisplacement,
5151 mat_grad_pos_at_pts_ptr));
5154 fe_rhs->getOpPtrVector().push_back(
5155 new HookeElement::OpCalculateStiffnessScaledByDensityField(
5156 "SPATIAL_POSITION",
"SPATIAL_POSITION", block_sets_ptr,
5158 fe_rhs->getOpPtrVector().push_back(
new HookeElement::OpCalculateStrainAle(
5159 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5160 fe_rhs->getOpPtrVector().push_back(
new HookeElement::OpCalculateStress<1>(
5161 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5162 fe_rhs->getOpPtrVector().push_back(
new HookeElement::OpCalculateEnergy(
5163 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5165 fe_rhs->getOpPtrVector().push_back(
new HookeElement::OpCalculateEshelbyStress(
5166 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5168 fe_rhs->getOpPtrVector().push_back(
new HookeElement::OpAleRhs_dX(
5169 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5170 fe_rhs->getOpPtrVector().push_back(
new HookeElement::OpAleRhs_dx(
5171 "SPATIAL_POSITION",
"SPATIAL_POSITION", data_hooke_element_at_pts));
5174 *data_hooke_element_at_pts,
mwlsApprox->rhoAtGaussPts,
5180 "SPATIAL_POSITION", space_grad_pos_at_pts_ptr, fe_lhs->singularElement,
5183 fe_lhs->singularElement, fe_lhs->detS, fe_lhs->invSJac));
5185 if (test_with_mwls) {
5186 if (
mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5188 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_lhs,
mwlsApprox,
5191 fe_lhs->getOpPtrVector().push_back(
5193 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_lhs,
mwlsApprox));
5194 fe_lhs->getOpPtrVector().push_back(
5196 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr, fe_lhs,
mwlsApprox,
5201 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr,
mwlsApprox->rhoAtGaussPts,
5203 fe_lhs,
mwlsApprox->singularInitialDisplacement,
5204 mat_grad_pos_at_pts_ptr));
5207 fe_lhs->getOpPtrVector().push_back(
5208 new HookeElement::OpCalculateStiffnessScaledByDensityField(
5209 "SPATIAL_POSITION",
"SPATIAL_POSITION", block_sets_ptr,
5212 fe_lhs->getOpPtrVector().push_back(
new HookeElement::OpCalculateStrainAle(
5213 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5214 fe_lhs->getOpPtrVector().push_back(
new HookeElement::OpCalculateStress<1>(
5215 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5216 fe_lhs->getOpPtrVector().push_back(
new HookeElement::OpCalculateEnergy(
5217 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5218 fe_lhs->getOpPtrVector().push_back(
new HookeElement::OpCalculateEshelbyStress(
5219 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5221 fe_lhs->getOpPtrVector().push_back(
new HookeElement::OpAleLhs_dX_dX<1>(
5222 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5223 fe_lhs->getOpPtrVector().push_back(
new HookeElement::OpAleLhsPre_dX_dx<1>(
5224 "MESH_NODE_POSITIONS",
"SPATIAL_POSITION", data_hooke_element_at_pts));
5225 fe_lhs->getOpPtrVector().push_back(
new HookeElement::OpAleLhs_dX_dx(
5226 "MESH_NODE_POSITIONS",
"SPATIAL_POSITION", data_hooke_element_at_pts));
5227 fe_lhs->getOpPtrVector().push_back(
new HookeElement::OpAleLhs_dx_dx<1>(
5228 "SPATIAL_POSITION",
"SPATIAL_POSITION", data_hooke_element_at_pts));
5229 fe_lhs->getOpPtrVector().push_back(
new HookeElement::OpAleLhs_dx_dX<1>(
5230 "SPATIAL_POSITION",
"MESH_NODE_POSITIONS", data_hooke_element_at_pts));
5232 fe_lhs->getOpPtrVector().push_back(
5233 new HookeElement::OpAleLhsWithDensity_dX_dX(
5234 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS",
5235 data_hooke_element_at_pts,
mwlsApprox->rhoAtGaussPts,
5238 fe_lhs->getOpPtrVector().push_back(
5239 new HookeElement::OpAleLhsWithDensity_dx_dX(
5240 "SPATIAL_POSITION",
"MESH_NODE_POSITIONS", data_hooke_element_at_pts,
5244 boost::shared_ptr<MatrixDouble> mat_singular_disp_ptr =
nullptr;
5246 mat_singular_disp_ptr = boost::shared_ptr<MatrixDouble>(
5249 fe_lhs->getOpPtrVector().push_back(
5251 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS",
5252 data_hooke_element_at_pts,
mwlsApprox->rhoAtGaussPts,
5255 fe_lhs->getOpPtrVector().push_back(
5257 "SPATIAL_POSITION",
"MESH_NODE_POSITIONS",
5258 data_hooke_element_at_pts,
mwlsApprox->rhoAtGaussPts,
5262 fe_lhs->getOpPtrVector().push_back(
5264 *data_hooke_element_at_pts,
mwlsApprox->rhoAtGaussPts,
5269 fe_lhs->getOpPtrVector().push_back(
5271 *data_hooke_element_at_pts,
mwlsApprox->rhoAtGaussPts,
5285 boost::shared_ptr<CrackFrontElement> fe_rhs;
5286 boost::shared_ptr<CrackFrontElement> fe_lhs;
5299 "There is no such test defined");
5304 CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
5305 CHKERR DMSetType(dm,
"MOFEM");
5320 Vec front_f, tangent_front_f;
5321 CHKERR DMCreateGlobalVector(dm, &front_f);
5322 CHKERR VecDuplicate(front_f, &tangent_front_f);
5323 CHKERR VecSetOption(front_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
5324 CHKERR VecSetOption(tangent_front_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
5326 if (
smootherFe->smootherData.frontF != PETSC_NULL)
5328 if (
smootherFe->smootherData.tangentFrontF != PETSC_NULL)
5332 smootherFe->smootherData.tangentFrontF = tangent_front_f;
5349 VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
5353 boost::shared_ptr<FEMethod>
null;
5371 if (
m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
5402 if (
m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
5427 "There is no such test defined");
5432 "MESH_NODE_POSITIONS");
5438 CHKERR VecSetOption(
f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
5439 CHKERR VecSetOption(q, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
5443 CHKERR VecCopy(q, q_copy);
5451 CHKERR VecGhostUpdateBegin(
f, INSERT_VALUES, SCATTER_FORWARD);
5452 CHKERR VecGhostUpdateEnd(
f, INSERT_VALUES, SCATTER_FORWARD);
5454 PetscBool calculate_fd = PETSC_FALSE;
5456 &calculate_fd, PETSC_NULL);
5459 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
5460 CHKERR SNESSetDM(snes, dm);
5463 CHKERR SNESAppendOptionsPrefix(snes,
"test_mwls_");
5465 if (calculate_fd == PETSC_FALSE) {
5466 char testing_options[] =
5467 "-test_mwls_snes_test_jacobian "
5468 "-test_mwls_snes_test_jacobian_display "
5469 "-test_mwls_snes_no_convergence_test -test_mwls_snes_atol 0 "
5470 "-test_mwls_snes_rtol 0 -test_mwls_snes_max_it 1 -test_mwls_pc_type "
5472 CHKERR PetscOptionsInsertString(NULL, testing_options);
5474 char testing_options[] =
5475 "-test_mwls_snes_no_convergence_test -test_mwls_snes_atol 0 "
5476 "-test_mwls_snes_rtol 0 -test_mwls_snes_max_it 1 -test_mwls_pc_type "
5478 CHKERR PetscOptionsInsertString(NULL, testing_options);
5483 CHKERR SNESSetFromOptions(snes);
5492 CHKERR SNESSolve(snes, PETSC_NULL, q);
5494 if (calculate_fd == PETSC_TRUE) {
5496 CHKERR MatNorm(
m, NORM_INFINITY, &nrm_m0);
5498 char testing_options_fd[] =
"-test_mwls_snes_fd";
5499 CHKERR PetscOptionsInsertString(NULL, testing_options_fd);
5503 CHKERR SNESSetJacobian(snes, fd_m, fd_m,
SnesMat, snes_ctx);
5504 CHKERR SNESSetFromOptions(snes);
5508 CHKERR SNESSolve(snes, NULL, q_copy);
5510 CHKERR MatAXPY(
m, -1, fd_m, SUBSET_NONZERO_PATTERN);
5513 CHKERR MatNorm(
m, NORM_INFINITY, &nrm_m);
5514 PetscPrintf(PETSC_COMM_WORLD,
"Matrix norms %3.4e %3.4e %3.4e\n", nrm_m0,
5515 nrm_m, nrm_m / nrm_m0);
5518 const double tol = 1e-6;
5521 "Difference between hand-calculated tangent matrix and finite "
5522 "difference matrix is too big");
5526 CHKERR SNESDestroy(&snes);
5536 "Elastic element instance not created");
5539 materialFe = boost::shared_ptr<NonlinearElasticElement>(
5543 for (std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
5545 sit !=
elasticFe->setOfBlocks.end(); sit++) {
5546 materialFe->setOfBlocks[sit->first] = sit->second;
5547 materialFe->setOfBlocks[sit->first].forcesOnlyOnEntitiesRow =
5550 materialFe->commonData.spatialPositions =
"SPATIAL_POSITION";
5551 materialFe->commonData.meshPositions =
"MESH_NODE_POSITIONS";
5566 boost::shared_ptr<map<int, BlockData>> block_sets_ptr(
5568 boost::shared_ptr<HookeElement::DataAtIntegrationPts>
5569 data_hooke_element_at_pts(
new HookeElement::DataAtIntegrationPts());
5573 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr(
new MatrixDouble());
5575 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
5577 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
5578 boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
5580 mat_grad_pos_at_pts_ptr =
5584 mat_grad_pos_at_pts_ptr));
5587 mat_grad_pos_at_pts_ptr));
5589 boost::shared_ptr<MatrixDouble> space_grad_pos_at_pts_ptr;
5591 space_grad_pos_at_pts_ptr =
5595 "SPATIAL_POSITION", space_grad_pos_at_pts_ptr,
5599 "SPATIAL_POSITION", space_grad_pos_at_pts_ptr,
5611 "MESH_NODE_POSITIONS",
materialFe->commonData));
5617 for (map<int, NonlinearElasticElement::BlockData>::iterator sit =
5619 sit !=
materialFe->setOfBlocks.end(); sit++) {
5622 "SPATIAL_POSITION", sit->second,
materialFe->commonData,
5626 "MESH_NODE_POSITIONS", sit->second,
materialFe->commonData));
5633 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
5634 mat_grad_pos_at_pts_ptr = data_hooke_element_at_pts->HMat;
5638 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
5640 space_grad_pos_at_pts_ptr = data_hooke_element_at_pts->hMat;
5645 "mwlsApprox not allocated");
5648 if (
mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5665 new HookeElement::OpCalculateStiffnessScaledByDensityField(
5666 "SPATIAL_POSITION",
"SPATIAL_POSITION", block_sets_ptr,
5670 new HookeElement::OpCalculateStrainAle(
"MESH_NODE_POSITIONS",
5671 "MESH_NODE_POSITIONS",
5672 data_hooke_element_at_pts));
5674 new HookeElement::OpCalculateStress<1>(
"MESH_NODE_POSITIONS",
5675 "MESH_NODE_POSITIONS",
5676 data_hooke_element_at_pts));
5680 new HookeElement::OpCalculateHomogeneousStiffness<0>(
5681 "SPATIAL_POSITION",
"SPATIAL_POSITION", block_sets_ptr,
5682 data_hooke_element_at_pts));
5684 new HookeElement::OpCalculateStrainAle(
"MESH_NODE_POSITIONS",
5685 "MESH_NODE_POSITIONS",
5686 data_hooke_element_at_pts));
5688 new HookeElement::OpCalculateStress<0>(
"MESH_NODE_POSITIONS",
5689 "MESH_NODE_POSITIONS",
5690 data_hooke_element_at_pts));
5693 new HookeElement::OpCalculateEnergy(
"MESH_NODE_POSITIONS",
5694 "MESH_NODE_POSITIONS",
5695 data_hooke_element_at_pts));
5698 new HookeElement::OpCalculateEshelbyStress(
"MESH_NODE_POSITIONS",
5699 "MESH_NODE_POSITIONS",
5700 data_hooke_element_at_pts));
5706 feMaterialRhs->getOpPtrVector().push_back(
new HookeElement::OpAleRhs_dX(
5707 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS",
5708 data_hooke_element_at_pts));
5713 *data_hooke_element_at_pts,
mwlsApprox->rhoAtGaussPts,
5722 "mwlsApprox not allocated");
5729 if (
mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5747 mat_grad_pos_at_pts_ptr,
5751 new HookeElement::OpAnalyticalInternalAleStrain_dX<0>(
5752 "MESH_NODE_POSITIONS", data_hooke_element_at_pts,
5754 boost::shared_ptr<MatrixDouble>(
5760 #ifdef __ANALITICAL_TRACTION__
5764 boost::shared_ptr<NeumannForcesSurface::MyTriangleFE>(
5770 if (meshset_manager_ptr->
checkMeshset(
"ANALITICAL_TRACTION")) {
5773 &cubit_meshset_ptr);
5781 "MESH_NODE_POSITIONS", faces,
5789 #endif //__ANALITICAL_TRACTION__
5799 "MESH_NODE_POSITIONS",
materialFe->commonData));
5806 for (
auto sit =
materialFe->setOfBlocks.begin();
5807 sit !=
materialFe->setOfBlocks.end(); sit++) {
5810 "SPATIAL_POSITION", sit->second,
materialFe->commonData,
5814 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS", sit->second,
5820 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
5821 mat_grad_pos_at_pts_ptr = data_hooke_element_at_pts->HMat;
5824 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
5826 space_grad_pos_at_pts_ptr = data_hooke_element_at_pts->hMat;
5831 "mwlsApprox not allocated");
5834 if (
mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5851 new HookeElement::OpCalculateStiffnessScaledByDensityField(
5852 "SPATIAL_POSITION",
"SPATIAL_POSITION", block_sets_ptr,
5857 new HookeElement::OpCalculateStrainAle(
"MESH_NODE_POSITIONS",
5858 "MESH_NODE_POSITIONS",
5859 data_hooke_element_at_pts));
5862 new HookeElement::OpCalculateStress<1>(
"MESH_NODE_POSITIONS",
5863 "MESH_NODE_POSITIONS",
5864 data_hooke_element_at_pts));
5867 new HookeElement::OpCalculateEnergy(
"MESH_NODE_POSITIONS",
5868 "MESH_NODE_POSITIONS",
5869 data_hooke_element_at_pts));
5871 new HookeElement::OpCalculateEshelbyStress(
5872 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS",
5873 data_hooke_element_at_pts));
5879 new HookeElement::OpAleLhs_dX_dX<1>(
"MESH_NODE_POSITIONS",
5880 "MESH_NODE_POSITIONS",
5881 data_hooke_element_at_pts));
5884 new HookeElement::OpAleLhsWithDensity_dX_dX(
5885 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS",
5886 data_hooke_element_at_pts,
mwlsApprox->rhoAtGaussPts,
5889 boost::shared_ptr<MatrixDouble> mat_singular_disp_ptr =
nullptr;
5892 mat_singular_disp_ptr = boost::shared_ptr<MatrixDouble>(
5897 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS",
5898 data_hooke_element_at_pts,
mwlsApprox->rhoAtGaussPts,
5900 mat_singular_disp_ptr));
5905 *data_hooke_element_at_pts,
mwlsApprox->rhoAtGaussPts,
5912 new HookeElement::OpCalculateHomogeneousStiffness<0>(
5913 "SPATIAL_POSITION",
"SPATIAL_POSITION", block_sets_ptr,
5914 data_hooke_element_at_pts));
5917 new HookeElement::OpCalculateStrainAle(
"MESH_NODE_POSITIONS",
5918 "MESH_NODE_POSITIONS",
5919 data_hooke_element_at_pts));
5921 new HookeElement::OpCalculateStress<0>(
"MESH_NODE_POSITIONS",
5922 "MESH_NODE_POSITIONS",
5923 data_hooke_element_at_pts));
5926 new HookeElement::OpCalculateEnergy(
"MESH_NODE_POSITIONS",
5927 "MESH_NODE_POSITIONS",
5928 data_hooke_element_at_pts));
5930 new HookeElement::OpCalculateEshelbyStress(
5931 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS",
5932 data_hooke_element_at_pts));
5938 new HookeElement::OpAleLhs_dX_dX<0>(
"MESH_NODE_POSITIONS",
5939 "MESH_NODE_POSITIONS",
5940 data_hooke_element_at_pts));
5946 "mwlsApprox not allocated");
5949 if (
mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
5967 space_grad_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
mwlsApprox,
5976 const int verb,
const bool debug) {
5980 "Elastic element instance not created");
5993 for (std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
5995 sit !=
elasticFe->setOfBlocks.end(); sit++) {
5997 smootherFe->setOfBlocks[0].tEts.merge(sit->second.tEts);
6004 smootherFe->commonData.spatialPositions =
"MESH_NODE_POSITIONS";
6005 smootherFe->commonData.meshPositions =
"NONE";
6008 smootherFe->feRhs.meshPositionsFieldName =
"NONE";
6009 smootherFe->feLhs.meshPositionsFieldName =
"NONE";
6020 "MESH_NODE_POSITIONS",
smootherFe->commonData));
6021 map<int, NonlinearElasticElement::BlockData>::iterator sit =
6024 "MESH_NODE_POSITIONS",
smootherFe->setOfBlocks.at(0),
6027 "MESH_NODE_POSITIONS", sit->second,
smootherFe->commonData,
6033 "MESH_NODE_POSITIONS",
smootherFe->commonData));
6035 "MESH_NODE_POSITIONS",
smootherFe->setOfBlocks.at(0),
6038 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS",
6040 smootherFe->smootherData,
"LAMBDA_CRACKFRONT_AREA_TANGENT"));
6046 mField,
"LAMBDA_CRACKFRONT_AREA_TANGENT"));
6048 Range contact_faces;
6055 boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>(
6057 for (
auto id : ids) {
6058 if (
id != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
6063 "LAMBDA_SURFACE" + boost::lexical_cast<std::string>(
id),
6069 boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>(
6076 boost::lexical_cast<std::string>(
6077 getInterface<CPMeshCut>()->getCrackSurfaceId()),
6081 auto get_edges_block_set = [
this]() {
6082 return getInterface<CPMeshCut>()->getEdgesBlockSet();
6084 if (get_edges_block_set()) {
6086 boost::shared_ptr<EdgeSlidingConstrains>(
6101 Range fix_material_nodes;
6103 mField,
"MESH_NODE_POSITIONS", fix_material_nodes);
6106 for (
auto f : *fields_ptr) {
6107 if (
f->getName().compare(0, 6,
"LAMBDA") == 0 &&
6108 f->getName() !=
"LAMBDA_ARC_LENGTH" &&
6109 f->getName() !=
"LAMBDA_CONTACT" &&
6110 f->getName() !=
"LAMBDA_CLOSE_CRACK") {
6145 boost::shared_ptr<GriffithForceElement::MyTriangleFEConstrainsDelta>(
6147 mField,
"LAMBDA_CRACKFRONT_AREA"));
6154 boost::shared_ptr<GriffithForceElement::MyTriangleFEConstrains>(
6156 mField,
"LAMBDA_CRACKFRONT_AREA",
6160 boost::shared_ptr<GriffithForceElement::MyTriangleFEConstrains>(
6162 mField,
"LAMBDA_CRACKFRONT_AREA",
6184 "Elastic element instance not created");
6193 boost::shared_ptr<map<int, BlockData>> block_sets_ptr(
6202 "SPATIAL_POSITION",
elasticFe->commonData,
6208 "MESH_NODE_POSITIONS",
elasticFe->commonData));
6215 map<int, NonlinearElasticElement::BlockData>::iterator sit =
6217 for (; sit !=
elasticFe->setOfBlocks.end(); sit++) {
6221 "SPATIAL_POSITION", sit->second,
elasticFe->commonData,
6226 "SPATIAL_POSITION",
"SPATIAL_POSITION", sit->second,
6231 "SPATIAL_POSITION",
"MESH_NODE_POSITIONS", sit->second,
6235 boost::shared_ptr<HookeElement::DataAtIntegrationPts>
6236 data_hooke_element_at_pts(
new HookeElement::DataAtIntegrationPts());
6239 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
6241 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr(
new MatrixDouble());
6244 mat_pos_at_pts_ptr));
6249 "mwlsApprox not allocated");
6252 boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
6253 mat_grad_pos_at_pts_ptr =
6257 mat_grad_pos_at_pts_ptr));
6260 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
6264 if (
mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
6267 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6272 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6276 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6280 new HookeElement::OpCalculateStiffnessScaledByDensityField(
6281 "SPATIAL_POSITION",
"SPATIAL_POSITION", block_sets_ptr,
6288 new HookeElement::OpCalculateStrainAle(
"SPATIAL_POSITION",
6290 data_hooke_element_at_pts));
6292 new HookeElement::OpCalculateStress<1>(
"SPATIAL_POSITION",
6294 data_hooke_element_at_pts));
6301 new HookeElement::OpAleLhs_dx_dx<1>(
"SPATIAL_POSITION",
6303 data_hooke_element_at_pts));
6305 new HookeElement::OpAleLhs_dx_dX<1>(
"SPATIAL_POSITION",
6306 "MESH_NODE_POSITIONS",
6307 data_hooke_element_at_pts));
6310 new HookeElement::OpAleLhsWithDensity_dx_dX(
6311 "SPATIAL_POSITION",
"MESH_NODE_POSITIONS",
6312 data_hooke_element_at_pts,
mwlsApprox->rhoAtGaussPts,
6315 boost::shared_ptr<MatrixDouble> mat_singular_disp_ptr =
nullptr;
6318 mat_singular_disp_ptr = boost::shared_ptr<MatrixDouble>(
6323 "SPATIAL_POSITION",
"MESH_NODE_POSITIONS",
6324 data_hooke_element_at_pts,
mwlsApprox->rhoAtGaussPts,
6326 mat_singular_disp_ptr));
6332 new HookeElement::OpCalculateHomogeneousStiffness<0>(
6333 "SPATIAL_POSITION",
"SPATIAL_POSITION", block_sets_ptr,
6334 data_hooke_element_at_pts));
6340 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
6344 new HookeElement::OpCalculateStrainAle(
"SPATIAL_POSITION",
6346 data_hooke_element_at_pts));
6348 new HookeElement::OpCalculateStress<0>(
"SPATIAL_POSITION",
6350 data_hooke_element_at_pts));
6357 new HookeElement::OpAleLhs_dx_dx<0>(
"SPATIAL_POSITION",
6359 data_hooke_element_at_pts));
6361 new HookeElement::OpAleLhs_dx_dX<0>(
"SPATIAL_POSITION",
6362 "MESH_NODE_POSITIONS",
6363 data_hooke_element_at_pts));
6369 "Material element instance not created");
6378 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr(
new MatrixDouble());
6381 mat_pos_at_pts_ptr));
6382 boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
6384 mat_grad_pos_at_pts_ptr =
6388 mat_grad_pos_at_pts_ptr));
6390 boost::shared_ptr<MatrixDouble> space_grad_pos_at_pts_ptr;
6392 space_grad_pos_at_pts_ptr =
6396 "SPATIAL_POSITION", space_grad_pos_at_pts_ptr,
6410 "MESH_NODE_POSITIONS",
materialFe->commonData));
6415 for (map<int, NonlinearElasticElement::BlockData>::iterator sit =
6417 sit !=
materialFe->setOfBlocks.end(); sit++) {
6420 "SPATIAL_POSITION", sit->second,
materialFe->commonData,
6424 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS", sit->second,
6428 "MESH_NODE_POSITIONS",
"SPATIAL_POSITION", sit->second,
6433 boost::shared_ptr<HookeElement::DataAtIntegrationPts>
6434 data_hooke_element_at_pts(
new HookeElement::DataAtIntegrationPts());
6439 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
6440 mat_grad_pos_at_pts_ptr = data_hooke_element_at_pts->HMat;
6445 "mwlsApprox not allocated");
6450 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
6454 if (
mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
6457 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6462 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6466 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6471 new HookeElement::OpCalculateStiffnessScaledByDensityField(
6472 "SPATIAL_POSITION",
"SPATIAL_POSITION", block_sets_ptr,
6475 space_grad_pos_at_pts_ptr = data_hooke_element_at_pts->hMat;
6477 new HookeElement::OpCalculateStrainAle(
"MESH_NODE_POSITIONS",
6478 "MESH_NODE_POSITIONS",
6479 data_hooke_element_at_pts));
6481 new HookeElement::OpCalculateStress<1>(
"MESH_NODE_POSITIONS",
6482 "MESH_NODE_POSITIONS",
6483 data_hooke_element_at_pts));
6485 new HookeElement::OpCalculateEnergy(
"MESH_NODE_POSITIONS",
6486 "MESH_NODE_POSITIONS",
6487 data_hooke_element_at_pts));
6489 new HookeElement::OpCalculateEshelbyStress(
6490 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS",
6491 data_hooke_element_at_pts));
6497 new HookeElement::OpAleLhs_dX_dX<1>(
"MESH_NODE_POSITIONS",
6498 "MESH_NODE_POSITIONS",
6499 data_hooke_element_at_pts));
6501 new HookeElement::OpAleLhsPre_dX_dx<1>(
"MESH_NODE_POSITIONS",
6503 data_hooke_element_at_pts));
6505 new HookeElement::OpAleLhs_dX_dx(
"MESH_NODE_POSITIONS",
6507 data_hooke_element_at_pts));
6510 new HookeElement::OpAleLhsWithDensity_dX_dX(
6511 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS",
6512 data_hooke_element_at_pts,
mwlsApprox->rhoAtGaussPts,
6515 boost::shared_ptr<MatrixDouble> mat_singular_disp_ptr =
nullptr;
6518 mat_singular_disp_ptr = boost::shared_ptr<MatrixDouble>(
6523 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS",
6524 data_hooke_element_at_pts,
mwlsApprox->rhoAtGaussPts,
6526 mat_singular_disp_ptr));
6531 *data_hooke_element_at_pts,
mwlsApprox->rhoAtGaussPts,
6538 *data_hooke_element_at_pts,
mwlsApprox->rhoAtGaussPts,
6546 new HookeElement::OpCalculateHomogeneousStiffness<0>(
6547 "SPATIAL_POSITION",
"SPATIAL_POSITION", block_sets_ptr,
6548 data_hooke_element_at_pts));
6554 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
6557 space_grad_pos_at_pts_ptr = data_hooke_element_at_pts->hMat;
6559 new HookeElement::OpCalculateStrainAle(
"MESH_NODE_POSITIONS",
6560 "MESH_NODE_POSITIONS",
6561 data_hooke_element_at_pts));
6563 new HookeElement::OpCalculateStress<0>(
"MESH_NODE_POSITIONS",
6564 "MESH_NODE_POSITIONS",
6565 data_hooke_element_at_pts));
6567 new HookeElement::OpCalculateEnergy(
"MESH_NODE_POSITIONS",
6568 "MESH_NODE_POSITIONS",
6569 data_hooke_element_at_pts));
6571 new HookeElement::OpCalculateEshelbyStress(
6572 "MESH_NODE_POSITIONS",
"MESH_NODE_POSITIONS",
6573 data_hooke_element_at_pts));
6579 new HookeElement::OpAleLhs_dX_dX<0>(
"MESH_NODE_POSITIONS",
6580 "MESH_NODE_POSITIONS",
6581 data_hooke_element_at_pts));
6583 new HookeElement::OpAleLhsPre_dX_dx<0>(
"MESH_NODE_POSITIONS",
6585 data_hooke_element_at_pts));
6587 new HookeElement::OpAleLhs_dX_dx(
"MESH_NODE_POSITIONS",
6589 data_hooke_element_at_pts));
6596 "mwlsApprox not allocated");
6599 if (
mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
6602 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6607 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6611 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
6620 space_grad_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
mwlsApprox,
6624 space_grad_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
mwlsApprox,
6632 DM dm,
const bool fix_crack_front,
const int verb,
const bool debug) {
6633 boost::shared_ptr<FEMethod>
null;
6637 Vec front_f, tangent_front_f;
6638 CHKERR DMCreateGlobalVector(dm, &front_f);
6639 CHKERR VecDuplicate(front_f, &tangent_front_f);
6640 CHKERR VecSetOption(front_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
6641 CHKERR VecSetOption(tangent_front_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
6643 if (
smootherFe->smootherData.frontF != PETSC_NULL)
6645 if (
smootherFe->smootherData.tangentFrontF != PETSC_NULL)
6649 smootherFe->smootherData.tangentFrontF = tangent_front_f;
6651 CHKERR PetscObjectReference((PetscObject)front_f);
6652 CHKERR PetscObjectReference((PetscObject)tangent_front_f);
6662 CHKERR PetscObjectReference((PetscObject)front_f);
6663 CHKERR PetscObjectReference((PetscObject)tangent_front_f);
6664 CHKERR VecDestroy(&front_f);
6665 CHKERR VecDestroy(&tangent_front_f);
6672 VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
6684 cerr <<
"Node " << endl;
6685 FieldEntity_multiIndex::index<Ent_mi_tag>::type::iterator eit, hi_eit;
6686 eit = field_ents->get<
Ent_mi_tag>().lower_bound(*vit);
6687 hi_eit = field_ents->get<
Ent_mi_tag>().upper_bound(*vit);
6688 for (; eit != hi_eit; ++eit) {
6689 cerr << **eit << endl;
6707 if (
m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
6737 if (
m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
6760 CHKERR DMCreateGlobalVector(dm, &q);
6763 CHKERR VecSetOption(
f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
6765 CHKERR DMCreateMatrix(dm, &
m);
6767 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
6774 CHKERR SNESDestroy(&snes);
6776 MatView(
m, PETSC_VIEWER_STDOUT_WORLD);
6787 DM dm,
const bool fix_crack_front,
const int verb,
const bool debug) {
6788 boost::shared_ptr<FEMethod>
null;
6793 if (
smootherFe->smootherData.frontF != PETSC_NULL)
6795 if (
smootherFe->smootherData.tangentFrontF != PETSC_NULL)
6798 smootherFe->smootherData.frontF = PETSC_NULL;
6799 smootherFe->smootherData.tangentFrontF = PETSC_NULL;
6810 cerr <<
"Node " << endl;
6811 FieldEntity_multiIndex::index<Ent_mi_tag>::type::iterator eit, hi_eit;
6812 eit = field_ents->get<
Ent_mi_tag>().lower_bound(*vit);
6813 hi_eit = field_ents->get<
Ent_mi_tag>().upper_bound(*vit);
6814 for (; eit != hi_eit; ++eit) {
6815 cerr << **eit << endl;
6831 if (
m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
6852 if (
m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
6872 CHKERR DMCreateGlobalVector(dm, &q);
6875 CHKERR VecSetOption(
f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
6877 CHKERR DMCreateMatrix(dm, &
m);
6879 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
6886 CHKERR SNESDestroy(&snes);
6888 MatView(
m, PETSC_VIEWER_STDOUT_WORLD);
6899 const bool fix_small_g,
6902 Range fix_material_nodes;
6913 "Number of fixed nodes %d, number of fixed crack front nodes %d "
6915 fix_material_nodes.size(),
6922 DM dm, boost::shared_ptr<FEMethod> arc_method,
6923 boost::shared_ptr<ArcLengthCtx> arc_ctx,
6924 const std::vector<int> &surface_ids,
const int verb,
const bool debug) {
6925 boost::shared_ptr<FEMethod>
null;
6944 for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
6947 auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
6948 auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
6949 for (; oit != hi_oit; oit++) {
6950 if (boost::typeindex::type_id_runtime(*oit) ==
6951 boost::typeindex::type_id<NeumannForcesSurface::OpNeumannForce>()) {
6959 for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
6966 auto oit = fit->second->getLoopFeMatRhs().getOpPtrVector().begin();
6967 auto hi_oit = fit->second->getLoopFeMatRhs().getOpPtrVector().end();
6968 for (; oit != hi_oit; oit++) {
6969 if (boost::typeindex::type_id_runtime(*oit) ==
6970 boost::typeindex::type_id<
6974 .
F = arc_ctx->F_lambda;
6981 for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
6985 auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
6986 auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
6987 for (; oit != hi_oit; oit++) {
6988 if (boost::typeindex::type_id_runtime(*oit) ==
6989 boost::typeindex::type_id<
6994 if (boost::typeindex::type_id_runtime(*oit) ==
6995 boost::typeindex::type_id<
7004 for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
7011 auto oit = fit->second->getLoopFeMatRhs().getOpPtrVector().begin();
7012 auto hi_oit = fit->second->getLoopFeMatRhs().getOpPtrVector().end();
7013 for (; oit != hi_oit; oit++) {
7014 if (boost::typeindex::type_id_runtime(*oit) ==
7015 boost::typeindex::type_id<
7019 .
F = arc_ctx->F_lambda;
7025 for (boost::ptr_map<string, EdgeForce>::iterator fit =
edgeForces->begin();
7027 auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
7028 auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
7029 for (; oit != hi_oit; oit++) {
7030 if (boost::typeindex::type_id_runtime(*oit) ==
7031 boost::typeindex::type_id<EdgeForce::OpEdgeForce>()) {
7037 for (boost::ptr_map<string, NodalForce>::iterator fit =
nodalForces->begin();
7039 auto oit = fit->second->getLoopFe().getOpPtrVector().begin();
7040 auto hi_oit = fit->second->getLoopFe().getOpPtrVector().end();
7041 for (; oit != hi_oit; oit++) {
7042 if (boost::typeindex::type_id_runtime(*oit) ==
7043 boost::typeindex::type_id<NodalForce::OpNodalForce>()) {
7049 Vec front_f, tangent_front_f;
7050 CHKERR DMCreateGlobalVector(dm, &front_f);
7051 CHKERR VecDuplicate(front_f, &tangent_front_f);
7052 CHKERR VecSetOption(front_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
7053 CHKERR VecSetOption(tangent_front_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
7057 if (
smootherFe->smootherData.frontF != PETSC_NULL)
7059 if (
smootherFe->smootherData.tangentFrontF != PETSC_NULL)
7063 smootherFe->smootherData.tangentFrontF = tangent_front_f;
7082 VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
7110 if (
m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
7150 if (
m.first != getInterface<CPMeshCut>()->getCrackSurfaceId()) {
7172 for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
7176 dm, fit->first.c_str(),
7177 boost::shared_ptr<FEMethod>(
surfaceForces, &fit->second->getLoopFe()),
7183 for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
7187 dm, fit->first.c_str(),
7189 &fit->second->getLoopFeLhs()),
7193 dm, fit->first.c_str(),
7195 &fit->second->getLoopFeMatRhs()),
7198 dm, fit->first.c_str(),
7200 &fit->second->getLoopFeMatLhs()),
7207 for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
7211 dm, fit->first.c_str(),
7212 boost::shared_ptr<FEMethod>(
surfacePressure, &fit->second->getLoopFe()),
7217 for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
7221 dm, fit->first.c_str(),
7223 &fit->second->getLoopFeLhs()),
7226 dm, fit->first.c_str(),
7228 &fit->second->getLoopFeMatRhs()),
7231 dm, fit->first.c_str(),
7233 &fit->second->getLoopFeMatLhs()),
7250 PETSC_NULL, PETSC_NULL);
7257 for (boost::ptr_map<string, EdgeForce>::iterator fit =
edgeForces->begin();
7260 dm, fit->first.c_str(),
7261 boost::shared_ptr<FEMethod>(
edgeForces, &fit->second->getLoopFe()),
7264 for (boost::ptr_map<string, NodalForce>::iterator fit =
nodalForces->begin();
7267 dm, fit->first.c_str(),
7268 boost::shared_ptr<FEMethod>(
nodalForces, &fit->second->getLoopFe()),
7281 const std::string msg) {
7285 "Elastic element instance not created");
7287 feEnergy = boost::make_shared<CrackFrontElement>(
7291 feEnergy->meshPositionsFieldName =
"NONE";
7292 feRhs->addToRule = 0;
7296 feEnergy->getOpPtrVector().push_back(
7298 "SPATIAL_POSITION",
elasticFe->commonData,
7300 feEnergy->getOpPtrVector().push_back(
7302 "MESH_NODE_POSITIONS",
elasticFe->commonData));
7305 map<int, NonlinearElasticElement::BlockData>::iterator sit =
7307 for (; sit !=
elasticFe->setOfBlocks.end(); sit++) {
7308 feEnergy->getOpPtrVector().push_back(
7317 boost::shared_ptr<map<int, BlockData>> block_sets_ptr(
7320 boost::shared_ptr<HookeElement::DataAtIntegrationPts>
7321 data_hooke_element_at_pts(
new HookeElement::DataAtIntegrationPts());
7322 boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr =
7323 data_hooke_element_at_pts->HMat;
7324 boost::shared_ptr<MatrixDouble> space_grad_pos_at_pts_ptr =
7325 data_hooke_element_at_pts->hMat;
7327 feEnergy->getOpPtrVector().push_back(
7329 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
7334 "mwlsApprox not allocated");
7336 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr(
new MatrixDouble());
7338 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
7339 feEnergy->getOpPtrVector().push_back(
7341 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
7343 if (
mwlsApprox->getUseGlobalBaseAtMaterialReferenceConfiguration())
7344 feEnergy->getOpPtrVector().push_back(
7346 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
feEnergy,
7349 feEnergy->getOpPtrVector().push_back(
7351 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
feEnergy,
7353 feEnergy->getOpPtrVector().push_back(
7355 mat_pos_at_pts_ptr, mat_grad_pos_at_pts_ptr,
feEnergy,
7358 feEnergy->getOpPtrVector().push_back(
7359 new HookeElement::OpCalculateStiffnessScaledByDensityField(
7360 "SPATIAL_POSITION",
"SPATIAL_POSITION", block_sets_ptr,
7363 feEnergy->getOpPtrVector().push_back(
7364 new HookeElement::OpCalculateStrainAle(
"SPATIAL_POSITION",
7366 data_hooke_element_at_pts));
7368 feEnergy->getOpPtrVector().push_back(
7369 new HookeElement::OpCalculateStress<1>(
"SPATIAL_POSITION",
7371 data_hooke_element_at_pts));
7374 feEnergy->getOpPtrVector().push_back(
7375 new HookeElement::OpCalculateHomogeneousStiffness<0>(
7376 "SPATIAL_POSITION",
"SPATIAL_POSITION", block_sets_ptr,
7377 data_hooke_element_at_pts));
7381 feEnergy->getOpPtrVector().push_back(
7383 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat,
7385 feEnergy->getOpPtrVector().push_back(
7386 new HookeElement::OpCalculateStrainAle(
"SPATIAL_POSITION",
7388 data_hooke_element_at_pts));
7389 feEnergy->getOpPtrVector().push_back(
7390 new HookeElement::OpCalculateStress<0>(
"SPATIAL_POSITION",
7392 data_hooke_element_at_pts));
7397 feEnergy->getOpPtrVector().push_back(
new HookeElement::OpCalculateEnergy(
7398 "SPATIAL_POSITION",
"SPATIAL_POSITION", data_hooke_element_at_pts,
7402 feEnergy->snes_ctx = SnesMethod::CTX_SNESNONE;
7407 PetscPrintf(PETSC_COMM_WORLD,
"Elastic energy %8.8e\n", msg.c_str(),
7410 PetscPrintf(PETSC_COMM_WORLD,
"%s Elastic energy %8.8e\n", msg.c_str(),
7425 CHKERR VecSetOption(
f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
7427 CHKERR VecGhostUpdateBegin(
f, INSERT_VALUES, SCATTER_FORWARD);
7428 CHKERR VecGhostUpdateEnd(
f, INSERT_VALUES, SCATTER_FORWARD);
7433 #ifdef __ANALITICAL_TRACTION__
7439 #endif // __ANALITICAL_TRACTION__
7442 CHKERR VecGhostUpdateBegin(
f, ADD_VALUES, SCATTER_REVERSE);
7443 CHKERR VecGhostUpdateEnd(
f, ADD_VALUES, SCATTER_REVERSE);
7444 CHKERR VecGhostUpdateBegin(
f, INSERT_VALUES, SCATTER_FORWARD);
7445 CHKERR VecGhostUpdateEnd(
f, INSERT_VALUES, SCATTER_FORWARD);
7447 CHKERR VecView(
f, PETSC_VIEWER_STDOUT_WORLD);
7452 CHKERR VecGhostUpdateBegin(
f, INSERT_VALUES, SCATTER_FORWARD);
7453 CHKERR VecGhostUpdateEnd(
f, INSERT_VALUES, SCATTER_FORWARD);
7457 "MESH_NODE_POSITIONS");
7458 CHKERR VecGhostUpdateBegin(q, INSERT_VALUES, SCATTER_FORWARD);
7459 CHKERR VecGhostUpdateEnd(q, INSERT_VALUES, SCATTER_FORWARD);
7461 ent_method_material_position, 0,
7481 PetscFunctionReturn(0);
7484 CHKERR VecSetOption(
f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
7486 CHKERR VecGhostUpdateBegin(
f, INSERT_VALUES, SCATTER_FORWARD);
7487 CHKERR VecGhostUpdateEnd(
f, INSERT_VALUES, SCATTER_FORWARD);
7494 CHKERR VecGhostUpdateBegin(
f, ADD_VALUES, SCATTER_REVERSE);
7495 CHKERR VecGhostUpdateEnd(
f, ADD_VALUES, SCATTER_REVERSE);
7496 CHKERR VecGhostUpdateBegin(
f, INSERT_VALUES, SCATTER_FORWARD);
7497 CHKERR VecGhostUpdateEnd(
f, INSERT_VALUES, SCATTER_FORWARD);
7499 CHKERR VecView(
f, PETSC_VIEWER_STDOUT_WORLD);
7505 CHKERR VecGhostUpdateBegin(
f, INSERT_VALUES, SCATTER_FORWARD);
7506 CHKERR VecGhostUpdateEnd(
f, INSERT_VALUES, SCATTER_FORWARD);
7508 ent_method_material_smoothing, 0,
7528 PetscFunctionReturn(0);
7529 auto get_tag_data = [
this](std::string tag_name) {
7536 auto griffith_force = get_tag_data(
"GRIFFITH_FORCE");
7537 auto smoothing_force = get_tag_data(
"SMOOTHING_FORCE");
7539 auto get_node_tensor = [](std::vector<double> &data) {
7541 &data[0], &data[1], &data[2]);
7543 auto t_griffith_force = get_node_tensor(griffith_force);
7544 auto t_smoothing_force = get_node_tensor(smoothing_force);
7548 auto get_magnitude = [](
auto t) {
7550 return sqrt(
t(
i) *
t(
i));
7552 const double griffith_mgn = get_magnitude(t_griffith_force);
7553 const double smoothing_mg = get_magnitude(t_smoothing_force);
7555 PetscPrintf(
mField.
get_comm(),
"Node smoothing force fraction %ld %6.4e\n",
7558 ++t_smoothing_force;
7564 DM dm_surface, DM dm_project,
const std::vector<int> &ids,
const int verb,
7591 Range contact_faces;
7596 for (
auto id : ids) {
7597 bool crack_surface =
7598 (
id == getInterface<CPMeshCut>()->getCrackSurfaceId()) ?
true :
false;
7604 "LAMBDA_SURFACE" + boost::lexical_cast<std::string>(
id),
7605 "MESH_NODE_POSITIONS");
7606 std::string fe_name = (
id == getInterface<CPMeshCut>()->getCrackSurfaceId())
7611 &surface_constrain.
feLhs);
7614 auto get_edges_block_set = [
this]() {
7615 return getInterface<CPMeshCut>()->getEdgesBlockSet();
7617 if (get_edges_block_set()) {
7620 "MESH_NODE_POSITIONS");
7638 DM dm_front, DM dm_project,
const int verb,
const bool debug) {
7644 PetscPrintf(PETSC_COMM_WORLD,
"No DOFs at crack front\n");
7660 CHKERR MatShellSetOperation(Q, MATOP_MULT,
7670 constant_area.
feLhsPtr = boost::shared_ptr<ConstantArea::MyTriangleFE>(
7672 "LAMBDA_CRACKFRONT_AREA"));
7675 constant_area.
feLhsPtr->petscQ = Q;
7677 constant_area.
feLhsPtr->petscQ = PETSC_NULL;
7679 constant_area.
feLhsPtr->getOpPtrVector().push_back(
7682 constant_area.
feLhsPtr->getOpPtrVector().push_back(
7687 tangent_constrain.
feLhsPtr = boost::shared_ptr<ConstantArea::MyTriangleFE>(
7689 "LAMBDA_CRACKFRONT_AREA_TANGENT"));
7692 tangent_constrain.
feLhsPtr->petscQ = Q;
7694 tangent_constrain.
feLhsPtr->petscQ = PETSC_NULL;
7696 tangent_constrain.
feLhsPtr->getOpPtrVector().push_back(
7699 tangent_constrain.
feLhsPtr->getOpPtrVector().push_back(
7717 PetscViewerASCIIOpen(PETSC_COMM_WORLD,
"C.m", &viewer);
7718 PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
7720 PetscViewerPopFormat(viewer);
7721 PetscViewerDestroy(&viewer);
7732 DM dm,
Vec f_material_proj,
Vec f_griffith_proj,
Vec f_lambda,
7733 const double gc,
const int verb,
const bool debug) {
7737 PetscPrintf(PETSC_COMM_WORLD,
"No crack front projection data structure\n");
7747 CHKERR VecGetSize(f_material_proj, &
N);
7748 CHKERR VecGetLocalSize(f_material_proj, &
n);
7749 CHKERR VecGetSize(f_lambda, &
M);
7750 CHKERR VecGetLocalSize(f_lambda, &
m);
7753 CHKERR MatShellSetOperation(RT, MATOP_MULT,
7757 CHKERR MatMult(RT, f_material_proj, f_lambda);
7758 CHKERR VecScale(f_lambda, -1);
7759 CHKERR VecGhostUpdateBegin(f_lambda, INSERT_VALUES, SCATTER_FORWARD);
7760 CHKERR VecGhostUpdateEnd(f_lambda, INSERT_VALUES, SCATTER_FORWARD);
7765 CHKERR VecView(f_lambda, PETSC_VIEWER_STDOUT_WORLD);
7781 CHKERR VecGetArray(f_lambda, &a_lambda);
7788 for (; it != hi_it; ++it) {
7789 double val = a_lambda[it->get()->getPetscLocalDofIdx()];
7790 if (!std::isnormal(val)) {
7792 "Value vector is not a number val = %3.4f", val);
7794 if (it->get()->getName() ==
"LAMBDA_CRACKFRONT_AREA") {
7795 mapG1[it->get()->getEnt()] = val;
7796 }
else if (it->get()->getName() ==
"LAMBDA_CRACKFRONT_AREA_TANGENT") {
7797 mapG3[it->get()->getEnt()] = val;
7800 "Should not happen");
7803 CHKERR VecRestoreArray(f_lambda, &a_lambda);
7806 Vec f_mat_scat, f_griffith_scat;
7809 CHKERR VecDuplicate(f_mat_scat, &f_griffith_scat);
7811 INSERT_VALUES, SCATTER_FORWARD);
7813 INSERT_VALUES, SCATTER_FORWARD);
7815 f_griffith_scat, INSERT_VALUES, SCATTER_FORWARD);
7817 f_griffith_scat, INSERT_VALUES, SCATTER_FORWARD);
7818 CHKERR VecGhostUpdateBegin(f_mat_scat, INSERT_VALUES, SCATTER_FORWARD);
7819 CHKERR VecGhostUpdateEnd(f_mat_scat, INSERT_VALUES, SCATTER_FORWARD);
7820 CHKERR VecGhostUpdateBegin(f_griffith_scat, INSERT_VALUES, SCATTER_FORWARD);
7821 CHKERR VecGhostUpdateEnd(f_griffith_scat, INSERT_VALUES, SCATTER_FORWARD);
7822 double *a_material, *a_griffith;
7823 CHKERR VecGetArray(f_mat_scat, &a_material);
7824 CHKERR VecGetArray(f_griffith_scat, &a_griffith);
7831 for (; it != hi_it; it++) {
7832 if (it->get()->getName() ==
"MESH_NODE_POSITIONS") {
7836 a_material[it->get()->getPetscLocalDofIdx()];
7839 a_griffith[it->get()->getPetscLocalDofIdx()];
7842 "Should not happen");
7845 CHKERR VecRestoreArray(f_mat_scat, &a_material);
7846 CHKERR VecRestoreArray(f_griffith_scat, &a_griffith);
7847 CHKERR VecDestroy(&f_mat_scat);
7848 CHKERR VecDestroy(&f_griffith_scat);
7852 for (map<EntityHandle, double>::iterator mit =
mapG1.begin();
7853 mit !=
mapG1.end(); mit++) {
7855 double g1 = mit->second;
7857 double g3 =
mapG3[ent];
7862 double g2 = sqrt(fabs(
j *
j - g1 * g1 - g3 * g3));
7866 ss <<
"griffith force at ent ";
7867 ss << setw(5) << ent;
7870 ss <<
" " << setw(10) << setprecision(4) << coords[0];
7871 ss <<
" " << setw(10) << setprecision(4) << coords[1];
7872 ss <<
" " << setw(10) << setprecision(4) << coords[2];
7874 ss <<
"\t\tg1 " << scientific << setprecision(6) << g1;
7875 ss <<
" / " << scientific << setprecision(6) <<
j;
7876 ss <<
" ( " << scientific << setprecision(6) << g1 /
j <<
" )";
7878 ss <<
"\t\tg2 " << scientific << setprecision(6) << g2;
7879 ss <<
" ( " << scientific << setprecision(6) << g2 /
j <<
" )";
7881 ss <<
"\t\tg3 " << scientific << setprecision(6) << g3;
7882 ss <<
" ( " << scientific << setprecision(6) << g3 /
j <<
" )";
7884 ss <<
"\t relative error (gc-g)/gc ";
7885 ss << scientific << setprecision(6) << (gc - g1) / gc;
7886 ss <<
" / " << scientific << setprecision(6) << (gc -
j) / gc;
7889 MOFEM_LOG(
"CPWorld", Sev::inform) << ss.str();
7900 CHKERR VecDestroy(&vec_max);
7917 CHKERR VecView(
f, PETSC_VIEWER_STDOUT_WORLD);
7927 CHKERR VecView(x, PETSC_VIEWER_STDOUT_WORLD);
7932 CHKERR VecView(Cf, PETSC_VIEWER_STDOUT_WORLD);
7944 CHKERR MatShellSetOperation(Q, MATOP_MULT,
7952 "MATERIAL_FORCE_PROJECTED");
7953 CHKERR VecGhostUpdateBegin(f_proj, INSERT_VALUES, SCATTER_FORWARD);
7954 CHKERR VecGhostUpdateEnd(f_proj, INSERT_VALUES, SCATTER_FORWARD);
7956 ent_method_spatial_position, 0,
7970 CHKERR VecSetOption(f_griffith, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
7972 boost::scoped_ptr<GriffithForceElement> griffith_force_element(
7974 griffith_force_element->blockData[0].gc = gc;
7977 griffith_force_element->blockData[0]);
7978 gC = griffith_force_element->blockData[0].gc;
7980 griffith_force_element->feRhs.getOpPtrVector().push_back(
7983 griffith_force_element->commonData));
7984 griffith_force_element->feRhs.getOpPtrVector().push_back(
7986 griffith_force_element->blockData[0],
7987 griffith_force_element->commonData));
7988 griffith_force_element->feRhs.snes_f = f_griffith;
7989 CHKERR VecZeroEntries(f_griffith);
7990 CHKERR VecGhostUpdateBegin(f_griffith, INSERT_VALUES, SCATTER_FORWARD);
7991 CHKERR VecGhostUpdateEnd(f_griffith, INSERT_VALUES, SCATTER_FORWARD);
7993 &griffith_force_element->feRhs);
7994 CHKERR VecAssemblyBegin(f_griffith);
7995 CHKERR VecAssemblyEnd(f_griffith);
7996 CHKERR VecGhostUpdateBegin(f_griffith, ADD_VALUES, SCATTER_REVERSE);
7997 CHKERR VecGhostUpdateEnd(f_griffith, ADD_VALUES, SCATTER_REVERSE);
7998 CHKERR VecGhostUpdateBegin(f_griffith, INSERT_VALUES, SCATTER_FORWARD);
7999 CHKERR VecGhostUpdateEnd(f_griffith, INSERT_VALUES, SCATTER_FORWARD);
8003 CHKERR VecGhostUpdateBegin(f_griffith, INSERT_VALUES, SCATTER_FORWARD);
8004 CHKERR VecGhostUpdateEnd(f_griffith, INSERT_VALUES, SCATTER_FORWARD);
8012 Vec f_griffith_proj,
8024 double nrm2_griffith_force;
8025 CHKERR VecNorm(f_griffith, NORM_2, &nrm2_griffith_force);
8027 nrm2_griffith_force);
8031 CHKERR VecGetSize(f_griffith, &
M);
8032 CHKERR VecGetLocalSize(f_griffith, &
m);
8036 CHKERR MatShellSetOperation(Q, MATOP_MULT,
8038 CHKERR MatMult(Q, f_griffith, f_griffith_proj);
8041 CHKERR VecCopy(f_griffith, f_griffith_proj);
8047 CHKERR VecGhostUpdateBegin(f_griffith_proj, INSERT_VALUES, SCATTER_FORWARD);
8048 CHKERR VecGhostUpdateEnd(f_griffith_proj, INSERT_VALUES, SCATTER_FORWARD);
8050 mField, f_griffith_proj,
"GRIFFITH_FORCE_PROJECTED");
8052 ent_method_prj_griffith_force);
8054 double nrm2_griffith_force;
8055 CHKERR VecNorm(f_griffith_proj, NORM_2, &nrm2_griffith_force);
8057 nrm2_griffith_force);
8068 if (q != PETSC_NULL) {
8069 CHKERR VecGhostUpdateBegin(q, INSERT_VALUES, SCATTER_FORWARD);
8070 CHKERR VecGhostUpdateEnd(q, INSERT_VALUES, SCATTER_FORWARD);
8075 "SPATIAL_POSITION");
8077 ent_method_spatial_position, 0,
8084 std::ostringstream file_name;
8085 file_name <<
"out_crack_spatial_position.h5m";
8087 "PARALLEL=WRITE_PART", &meshset, 1);
8097 auto make_file_name = [step](
const std::string prefix,
8098 const std::string surfix =
".vtk") {
8099 return prefix + boost::lexical_cast<std::string>(step) + surfix;
8102 auto write_file = [
this, make_file_name](
const std::string prefix,
8107 std::string file_name = make_file_name(prefix);
8121 if (!crack_front_edge.empty())
8122 CHKERR write_file(
"out_crack_front_", crack_front_edge);
8126 double edge_node_coords[6];
8129 &edge_node_coords[2]),
8131 &edge_node_coords[5])};
8135 double sum_crack_front_length = 0;
8136 for (
auto edge : crack_front_edge) {
8141 t_node_edge[0](
i) -= t_node_edge[1](
i);
8142 double l = sqrt(t_node_edge[0](
i) * t_node_edge[0](
i));
8143 sum_crack_front_length +=
l;
8146 sum_crack_front_length);
8148 double rate_cf = 0.;
8151 PetscPrintf(
mField.
get_comm(),
" | Change rate = %6.4e\n", rate_cf);
8161 CHKERR write_file(
"out_crack_surface_", crack_surface);
8171 CHKERR Tools::getTriNormal(coords, &t_normal(0));
8173 return sqrt(t_normal(
i) * t_normal(
i)) * 0.5;
8176 double sum_crack_area = 0;
8177 for (
auto tri : crack_surface) {
8178 sum_crack_area += get_tri_area(tri);
8180 PetscPrintf(
mField.
get_comm(),
"Crack surface area = %6.4e", sum_crack_area);
8182 double rate_cs = 0.;
8185 PetscPrintf(
mField.
get_comm(),
" | Change rate = %6.4e\n", rate_cs);
8195 CHKERR write_file(
"out_volume_", bit_tets);
8199 CHKERR skin.find_skin(0, bit_tets,
false, tets_skin);
8200 CHKERR write_file(
"out_skin_volume_", tets_skin);
8210 CHKERR VecPointwiseMult(
f,
f, arc_snes_ctx->getVecDiagM());
8220 CHKERR MatDiagonalScale(B, arc_snes_ctx->getVecDiagM(), PETSC_NULL);
8225 Vec f,
bool snes_set_up,
8232 auto check_vec_size = [&](
auto q) {
8235 CHKERR VecGetSize(q, &q_size);
8238 "MoFEM vector size incomatibile %d != d", q_size,
8243 CHKERR check_vec_size(q);
8246 decltype(
arcCtx) arc_ctx;
8247 if (problem_ptr->
getName() ==
"ELASTIC")
8249 if (problem_ptr->
getName() ==
"EIGEN_ELASTIC")
8254 "Ctx for arc-length is not created");
8261 arc_snes_ctx->setVecDiagM(diag_m);
8262 CHKERR VecSet(arc_snes_ctx->getVecDiagM(), 1.);
8265 boost::make_shared<ArcLengthMatShell>(
m, arc_ctx, problem_ptr->
getName());
8267 auto create_shell_matrix = [&]() {
8272 CHKERR MatGetLocalSize(
m, &mm, &nn);
8274 static_cast<void *
>(arc_mat_ctx.get()), shellM);
8275 CHKERR MatShellSetOperation(*shellM, MATOP_MULT,
8281 auto get_pc = [&]() {
8283 CHKERR PCAppendOptionsPrefix(pc,
"elastic_");
8284 CHKERR PCSetFromOptions(pc);
8288 boost::shared_ptr<PCArcLengthCtx> pc_arc_length_ctx;
8290 auto get_pc_arc = [&](
auto pc_mg) {
8292 boost::make_shared<PCArcLengthCtx>(pc_mg, *shellM,
m, arc_ctx);
8294 CHKERR PCSetType(pc_arc, PCSHELL);
8295 CHKERR PCShellSetContext(pc_arc, pc_arc_length_ctx.get());
8301 auto setup_snes = [&](
auto pc_mg,
auto pc_arc) {
8305 CHKERR SNESSetDM(snes, dm);
8308 CHKERR SNESAppendOptionsPrefix(snes,
"elastic_");
8309 CHKERR SNESSetFromOptions(snes);
8310 CHKERR SNESGetKSP(snes, &ksp);
8311 CHKERR KSPSetPC(ksp, pc_arc);
8313 #if PETSC_VERSION_GE(3, 8, 0)
8314 CHKERR PCFactorSetMatSolverType(pc_mg, MATSOLVERMUMPS);
8316 CHKERR PCFactorSetMatSolverPackage(pc_mg, MATSOLVERMUMPS);
8323 auto calculate_diagonal_scale = [&]() {
8329 const double *f_array;
8330 CHKERR VecGetArrayRead(
f, &f_array);
8335 constexpr
size_t nb_fields = 2;
8336 std::array<std::string, nb_fields> fields = {
"SPATIAL_POSITION",
8337 "LAMBDA_CLOSE_CRACK"};
8338 std::array<double, nb_fields> lnorms = {0, 0};
8339 std::array<double, nb_fields> norms = {0, 0};
8340 std::array<double, nb_fields> loc_nb_dofs = {0, 0};
8341 std::array<double, nb_fields> glob_nb_dofs = {0, 0};
8345 for (
size_t f = 0;
f != nb; ++
f) {
8348 rows_index.lower_bound(FieldEntity::getLoBitNumberUId(bit_number));
8350 rows_index.upper_bound(FieldEntity::getHiBitNumberUId(bit_number));
8351 for (
auto lo = lo_uid; lo != hi_uid; ++lo) {
8352 const int local_idx = (*lo)->getPetscLocalDofIdx();
8353 if (local_idx >= 0 && local_idx < nb_local) {
8354 const double val = f_array[local_idx];
8355 lnorms[
f] += val * val;
8361 CHKERR VecRestoreArrayRead(
f, &f_array);
8362 CHKERR MPIU_Allreduce(lnorms.data(), norms.data(), nb, MPIU_REAL, MPIU_SUM,
8363 PetscObjectComm((PetscObject)dm));
8364 CHKERR MPIU_Allreduce(loc_nb_dofs.data(), glob_nb_dofs.data(), nb,
8365 MPIU_REAL, MPIU_SUM,
8366 PetscObjectComm((PetscObject)dm));
8368 for (
auto &
v : norms)
8371 for (
size_t f = 0;
f != nb; ++
f)
8373 "Scaling diaginal for field [ %s ] by %9.8e nb. of dofs %d",
8374 fields[
f].c_str(), norms[
f],
8375 static_cast<double>(glob_nb_dofs[
f]));
8377 CHKERR VecSet(arc_snes_ctx->getVecDiagM(), 1.);
8378 double *diag_m_array;
8379 CHKERR VecGetArray(arc_snes_ctx->getVecDiagM(), &diag_m_array);
8381 for (
size_t f = 0;
f != nb; ++
f) {
8384 const auto lo_uid = FieldEntity::getLoBitNumberUId(bit_number);
8385 const auto hi_uid = FieldEntity::getHiBitNumberUId(bit_number);
8386 for (
auto lo = rows_index.lower_bound(lo_uid);
8387 lo != rows_index.upper_bound(hi_uid); ++lo) {
8388 const int local_idx = (*lo)->getPetscLocalDofIdx();
8389 if (local_idx >= 0 && local_idx < nb_local)
8390 diag_m_array[local_idx] = 1. / norms[
f];
8395 CHKERR VecGetArray(arc_snes_ctx->getVecDiagM(), &diag_m_array);
8401 CHKERR create_shell_matrix();
8402 auto pc_base = get_pc();
8403 auto pc_arc = get_pc_arc(pc_base);
8404 CHKERR setup_snes(pc_base, pc_arc);
8409 CHKERR VecGhostUpdateBegin(
f, INSERT_VALUES, SCATTER_FORWARD);
8410 CHKERR VecGhostUpdateEnd(
f, INSERT_VALUES, SCATTER_FORWARD);
8418 CHKERR calculate_diagonal_scale();
8420 auto create_section = [&]() {
8422 PetscSection section;
8423 CHKERR DMGetDefaultSection(dm, §ion);
8425 CHKERR PetscSectionGetNumFields(section, &num_fields);
8426 for (
int f = 0;
f != num_fields;
f++) {
8435 auto set_monitor = [&]() {
8437 PetscViewerAndFormat *vf;
8438 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
8439 PETSC_VIEWER_DEFAULT, &vf);
8451 CHKERR SNESSolve(snes, PETSC_NULL, q);
8452 CHKERR VecGhostUpdateBegin(q, INSERT_VALUES, SCATTER_FORWARD);
8453 CHKERR VecGhostUpdateEnd(q, INSERT_VALUES, SCATTER_FORWARD);
8454 CHKERR SNESMonitorCancel(snes);
8457 CHKERR VecNorm(q, NORM_2, &fnorm);
8458 MOFEM_LOG_C(
"CPWorld", Sev::verbose,
"solution fnorm = %9.8e", fnorm);
8460 if (problem_ptr->
getName() ==
"EIGEN_ELASTIC") {
8462 "EIGEN_ELASTIC",
"SPATIAL_POSITION",
"EIGEN_SPATIAL_POSITIONS",
ROW, q,
8463 INSERT_VALUES, SCATTER_REVERSE);
8474 CHKERR VecPointwiseMult(
f,
f, arc_snes_ctx->getVecDiagM());
8486 CHKERR MatDiagonalScale(B, arc_snes_ctx->getVecDiagM(), PETSC_NULL);
8487 CHKERR VecPointwiseMult(arc_snes_ctx->getArcPtr()->F_lambda,
8488 arc_snes_ctx->getArcPtr()->F_lambda,
8489 arc_snes_ctx->getVecDiagM());
8495 SNES snes, Mat
m,
Vec q,
8500 Mat m_elastic = PETSC_NULL;
8511 arc_snes_ctx->setVecDiagM(diag_m);
8512 CHKERR VecSet(arc_snes_ctx->getVecDiagM(), 1.);
8519 mwlsApprox->F_lambda = arc_snes_ctx->getArcPtr()->F_lambda;
8525 for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
8533 for (boost::ptr_map<string, NeumannForcesSurface>::iterator fit =
8540 auto arc_mat_ctx = boost::make_shared<ArcLengthMatShell>(
8541 m, arc_snes_ctx->getArcPtr(), problem_ptr->
getName());
8543 auto create_shell_matrix = [&]() {
8548 CHKERR MatGetLocalSize(
m, &mm, &nn);
8550 static_cast<void *
>(arc_mat_ctx.get()), &shell_m);
8551 CHKERR MatShellSetOperation(shell_m, MATOP_MULT,
8556 auto set_monitor = [&]() {
8558 PetscViewerAndFormat *vf;
8559 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
8560 PETSC_VIEWER_DEFAULT, &vf);
8569 auto zero_matrices_and_vectors = [&]() {
8573 CHKERR VecGhostUpdateBegin(
f, INSERT_VALUES, SCATTER_FORWARD);
8574 CHKERR VecGhostUpdateEnd(
f, INSERT_VALUES, SCATTER_FORWARD);
8578 auto create_section = [&]() {
8580 PetscSection section;
8581 CHKERR DMGetDefaultSection(dm, §ion);
8583 CHKERR PetscSectionGetNumFields(section, &num_fields);
8584 for (
int f = 0;
f != num_fields;
f++) {
8593 boost::shared_ptr<PCArcLengthCtx> pc_arc_length_ctx;
8595 auto create_pc_arc = [&](
auto pc_fs) {
8596 pc_arc_length_ctx = boost::make_shared<PCArcLengthCtx>(
8597 pc_fs, shell_m,
m, arc_snes_ctx->getArcPtr());
8599 CHKERR PCSetType(pc_arc, PCSHELL);
8600 CHKERR PCShellSetContext(pc_arc, pc_arc_length_ctx.get());
8606 auto set_up_snes = [&](
auto pc_arc,
auto pc_fs,
auto is_pcfs) {
8611 CHKERR SNESSetDM(snes, dm);
8614 CHKERR SNESSetFromOptions(snes);
8617 CHKERR SNESGetKSP(snes, &ksp);
8618 CHKERR KSPSetPC(ksp, pc_arc);
8622 CHKERR PCFieldSplitGetSubKSP(pc_fs, &
n, &sub_ksp);
8624 CHKERR KSPGetPC(sub_ksp[0], &pc_0);
8625 CHKERR KSPGetPC(sub_ksp[1], &pc_1);
8626 #if PETSC_VERSION_GE(3, 8, 0)
8627 CHKERR PCFactorSetMatSolverType(pc_0, MATSOLVERMUMPS);
8628 CHKERR PCFactorSetMatSolverType(pc_1, MATSOLVERMUMPS);
8630 CHKERR PCFactorSetMatSolverPackage(pc_0, MATSOLVERMUMPS);
8631 CHKERR PCFactorSetMatSolverPackage(pc_1, MATSOLVERMUMPS);
8635 #if PETSC_VERSION_GE(3, 8, 0)
8636 CHKERR PCFactorSetMatSolverType(pc_fs, MATSOLVERMUMPS);
8638 CHKERR PCFactorSetMatSolverPackage(pc_fs, MATSOLVERMUMPS);
8645 auto create_pc_fs = [&](
auto &is_pcfs) {
8647 CHKERR PCAppendOptionsPrefix(pc_fs,
"propagation_");
8648 CHKERR PCSetFromOptions(pc_fs);
8649 PetscObjectTypeCompare((PetscObject)pc_fs, PCFIELDSPLIT, &is_pcfs);
8653 boost::shared_ptr<ComposedProblemsData> cmp_data_ptr =
8655 for (
int ff = 0; ff !=
static_cast<int>(cmp_data_ptr->rowIs.size());
8657 CHKERR PCFieldSplitSetIS(pc_fs, NULL, cmp_data_ptr->rowIs[ff]);
8659 CHKERR PCSetOperators(pc_fs,
m,
m);
8665 auto calculate_diagonal_scale = [&](
bool debug =
false) {
8675 auto f_lambda = arc_snes_ctx->getArcPtr()->F_lambda;
8682 CHKERR VecSetOption(f_griffith, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
8684 auto assemble_mat = [&](DM dm, std::string name,
8685 boost::shared_ptr<FEMethod> fe, Mat
m) {
8687 fe->snes_ctx = SnesMethod::CTX_SNESSETJACOBIAN;
8691 fe->snes_ctx = SnesMethod::CTX_SNESNONE;
8695 auto assemble_vec = [&](DM dm, std::string name,
8696 boost::shared_ptr<FEMethod> fe,
Vec f) {
8698 fe->snes_ctx = SnesMethod::CTX_SNESSETFUNCTION;
8702 fe->snes_ctx = SnesMethod::CTX_SNESNONE;
8706 auto calculate_mat_and_vec = [&](
double set_arc_alpha,
8707 double set_smoothing_alpha,
double set_gc,
8710 bool stabilised =
smootherFe->smootherData.sTabilised;
8711 double arc_alpha = arc_snes_ctx->getArcPtr()->alpha;
8719 arc_snes_ctx->getArcPtr()->alpha = set_arc_alpha;
8725 CHKERR VecZeroEntries(f_griffith);
8728 CHKERR VecGhostUpdateBegin(f_griffith, ADD_VALUES, SCATTER_REVERSE);
8729 CHKERR VecGhostUpdateEnd(f_griffith, ADD_VALUES, SCATTER_REVERSE);
8730 CHKERR VecAssemblyBegin(f_griffith);
8731 CHKERR VecAssemblyEnd(f_griffith);
8735 CHKERR MatAssemblyBegin(
m, MAT_FINAL_ASSEMBLY);
8736 CHKERR MatAssemblyEnd(
m, MAT_FINAL_ASSEMBLY);
8737 CHKERR MatGetDiagonal(
m, diag_griffith);
8742 CHKERR MatAssemblyBegin(
m, MAT_FINAL_ASSEMBLY);
8743 CHKERR MatAssemblyEnd(
m, MAT_FINAL_ASSEMBLY);
8744 CHKERR MatGetDiagonal(
m, diag_smoothing);
8746 smootherFe->smootherData.sTabilised = stabilised;
8747 arc_snes_ctx->getArcPtr()->alpha = arc_alpha;
8755 auto get_norm_from_vec = [&](
Vec vec, std::string field,
Range *ents,
8759 const double *v_array;
8760 CHKERR VecGetArrayRead(vec, &v_array);
8763 const auto lo_uid = FieldEntity::getLoBitNumberUId(bit_number);
8764 const auto hi_uid = FieldEntity::getHiBitNumberUId(bit_number);
8766 for (
auto lo = rows_index.lower_bound(lo_uid);
8767 lo != rows_index.upper_bound(hi_uid); ++lo) {
8769 const int local_idx = (*lo)->getPetscLocalDofIdx();
8770 if (local_idx >= 0 && local_idx < nb_local) {
8772 bool no_skip =
true;
8774 bool contain = ents->find((*lo)->getEnt()) != ents->end();
8785 const double val = v_array[local_idx];
8786 ret_val += val * val;
8791 CHKERR VecRestoreArrayRead(vec, &v_array);
8792 return std::pair<double, double>(ret_val, size);
8795 auto set_norm_from_vec = [&](
Vec vec, std::string field,
const double norm,
8796 Range *ents,
bool not_contain) {
8799 if (!std::isnormal(norm))
8801 "Wrong scaling value");
8804 CHKERR VecGetArray(vec, &v_array);
8807 const auto lo_uid = FieldEntity::getLoBitNumberUId(bit_number);
8808 const auto hi_uid = FieldEntity::getHiBitNumberUId(bit_number);
8810 for (
auto lo = rows_index.lower_bound(lo_uid);
8811 lo != rows_index.upper_bound(hi_uid); ++lo) {
8813 const int local_idx = (*lo)->getPetscLocalDofIdx();
8814 if (local_idx >= 0 && local_idx < nb_local) {
8816 bool no_skip =
true;
8818 bool contain = ents->find((*lo)->getEnt()) != ents->end();
8829 v_array[local_idx] /= norm;
8832 CHKERR VecRestoreArray(vec, &v_array);
8836 auto calulate_norms = [&](
auto &lnorms,
auto &lsizes,
auto &norms,
8839 const double lambda = arc_snes_ctx->getArcPtr()->getFieldData();
8840 std::pair<double &, double &>(lnorms[0], lsizes[0]) =
8841 get_norm_from_vec(
f,
"LAMBDA_ARC_LENGTH",
nullptr,
false);
8842 std::pair<double &, double &>(lnorms[1], lsizes[1]) =
8843 get_norm_from_vec(f_lambda,
"SPATIAL_POSITION",
nullptr,
false);
8844 std::pair<double &, double &>(lnorms[2], lsizes[2]) = get_norm_from_vec(
8846 std::pair<double &, double &>(lnorms[3], lsizes[3]) = get_norm_from_vec(
8848 std::pair<double &, double &>(lnorms[4], lsizes[4]) = get_norm_from_vec(
8851 CHKERR MPIU_Allreduce(lnorms.data(), norms.data(), lnorms.size(),
8852 MPIU_REAL, MPIU_SUM,
8853 PetscObjectComm((PetscObject)dm));
8854 CHKERR MPIU_Allreduce(lsizes.data(), sizes.data(), lnorms.size(),
8855 MPIU_REAL, MPIU_SUM,
8856 PetscObjectComm((PetscObject)dm));
8858 for (
auto &
v : norms)
8862 norms[3] /= norms[1];
8863 norms[4] /= norms[1];
8864 norms[3] /= sizes[3];
8865 norms[4] /= sizes[4];
8870 auto test_scale = [&](
auto &fields,
auto &lnorms,
auto &lsizes,
auto &norms,
8873 CHKERR calculate_mat_and_vec(arc_snes_ctx->getArcPtr()->alpha,
8878 for (
auto &
v : lnorms)
8880 for (
auto &
v : lsizes)
8882 for (
auto &
v : norms)
8884 for (
auto &
v : sizes)
8887 CHKERR VecPointwiseMult(
f,
f, arc_snes_ctx->getVecDiagM());
8888 CHKERR VecPointwiseMult(f_lambda, f_lambda, arc_snes_ctx->getVecDiagM());
8889 CHKERR VecPointwiseMult(f_griffith, f_griffith,
8890 arc_snes_ctx->getVecDiagM());
8891 CHKERR VecPointwiseMult(diag_griffith, diag_griffith,
8892 arc_snes_ctx->getVecDiagM());
8893 CHKERR VecPointwiseMult(diag_smoothing, diag_smoothing,
8894 arc_snes_ctx->getVecDiagM());
8896 CHKERR calulate_norms(lnorms, lsizes, norms, sizes);
8901 for (
size_t f = 0;
f != fields.size(); ++
f)
8903 "Norms for field [ %s ] = %9.8e (%6.0f)", fields[
f].c_str(),
8904 norms[
f], sizes[
f]);
8908 CHKERR calculate_mat_and_vec(1, 1,
gC, 1);
8910 constexpr
size_t nb_fields = 5;
8911 std::array<std::string, nb_fields> fields = {
8912 "LAMBDA_ARC_LENGTH",
"FLAMBDA",
"GRIFFITH_FORCE",
"GRIFFITH_FORCE_LHS",
8913 "SMOOTHING_CONSTRAIN_LHS"};
8914 std::array<double, nb_fields> lnorms = {0, 0, 0, 0, 0};
8915 std::array<double, nb_fields> norms = {0, 0, 0, 0, 0};
8916 std::array<double, nb_fields> lsizes = {0, 0, 0, 0, 0};
8917 std::array<double, nb_fields> sizes = {0, 0, 0, 0, 0};
8919 CHKERR calulate_norms(lnorms, lsizes, norms, sizes);
8921 constexpr
size_t nb_fields_contact = 2;
8922 std::array<std::string, nb_fields> fields_contact = {
"LAMBDA_CONTACT"};
8923 std::array<double, nb_fields> lnorms_contact = {0, 0};
8924 std::array<double, nb_fields> norms_contact = {0, 0};
8925 std::array<double, nb_fields> lsizes_contact = {0, 0};
8926 std::array<double, nb_fields> sizes_contact = {0, 0};
8930 for (
size_t f = 0;
f != nb_fields; ++
f)
8932 "Norms for field [ %s ] = %9.8e (%6.0f)", fields[
f].c_str(),
8933 norms[
f], sizes[
f]);
8935 arc_snes_ctx->getArcPtr()->alpha =
arcAlpha / norms[0];
8936 const double scaled_smoother_alpha =
smootherAlpha * norms[3] / norms[4];
8940 CHKERR set_norm_from_vec(arc_snes_ctx->getVecDiagM(),
"SPATIAL_POSITION",
8941 norms[1],
nullptr,
false);
8942 CHKERR set_norm_from_vec(arc_snes_ctx->getVecDiagM(),
"MESH_NODE_POSITIONS",
8943 norms[1],
nullptr,
false);
8947 MOFEM_LOG_C(
"CPWorld", Sev::inform,
"Set scaled arc length alpha = %4.3e",
8948 arc_snes_ctx->getArcPtr()->alpha);
8949 MOFEM_LOG_C(
"CPWorld", Sev::inform,
"Set scaled smoothing alpha = %4.3e",
8951 MOFEM_LOG_C(
"CPWorld", Sev::inform,
"Set scaled griffith E = %4.3e",
8955 CHKERR test_scale(fields, lnorms, lsizes, norms, sizes);
8962 PetscBool is_pcfs = PETSC_FALSE;
8964 CHKERR create_shell_matrix();
8965 auto pc_fs = create_pc_fs(is_pcfs);
8966 auto pc_arc = create_pc_arc(pc_fs);
8967 CHKERR set_up_snes(pc_arc, pc_fs, is_pcfs);
8971 CHKERR zero_matrices_and_vectors();
8980 CHKERR calculate_diagonal_scale(
false);
8985 <<
"Resest MWLS approximation coefficients. Coefficients will be "
8986 "recalulated for current material positions.";
8993 CHKERR SNESSolve(snes, PETSC_NULL, q);
8994 CHKERR VecGhostUpdateBegin(q, INSERT_VALUES, SCATTER_FORWARD);
8995 CHKERR VecGhostUpdateEnd(q, INSERT_VALUES, SCATTER_FORWARD);
8996 CHKERR SNESMonitorCancel(snes);
8999 CHKERR VecNorm(q, NORM_2, &fnorm);
9002 CHKERR MatDestroy(&shell_m);
9004 CHKERR MatDestroy(&m_elastic);
9011 const std::string fe_name,
9012 const bool approx_internal_stress) {
9016 "Pointer to elasticFe is NULL");
9019 auto post_proc = boost::make_shared<
9028 post_proc->meshPositionsFieldName =
"NONE";
9030 CHKERR post_proc->generateReferenceElementMesh();
9031 CHKERR post_proc->addFieldValuesPostProc(
"MESH_NODE_POSITIONS");
9032 CHKERR post_proc->addFieldValuesGradientPostProc(
"MESH_NODE_POSITIONS");
9035 boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr(
new MatrixDouble());
9037 "MESH_NODE_POSITIONS", mat_pos_at_pts_ptr));
9038 boost::shared_ptr<MatrixDouble> mat_grad_pos_at_pts_ptr;
9041 post_proc->getOpPtrVector().push_back(
9044 post_proc->singularElement, post_proc->singularDisp,
9045 post_proc->postProcMesh, post_proc->mapGaussPts,
H));
9046 post_proc->getOpPtrVector().push_back(
9048 "MESH_NODE_POSITIONS",
elasticFe->commonData));
9050 post_proc->singularElement, post_proc->detS, post_proc->invSJac));
9051 CHKERR post_proc->addFieldValuesPostProc(
"SPATIAL_POSITION");
9052 CHKERR post_proc->addFieldValuesGradientPostProc(
"SPATIAL_POSITION");
9053 std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
9055 for (; sit !=
elasticFe->setOfBlocks.end(); sit++) {
9057 post_proc->postProcMesh, post_proc->mapGaussPts,
"SPATIAL_POSITION",
9058 sit->second, post_proc->commonData,
9064 boost::shared_ptr<HookeElement::DataAtIntegrationPts>
9065 data_hooke_element_at_pts(
new HookeElement::DataAtIntegrationPts());
9079 struct OpGetFieldGradientValuesOnSkinWithSingular
9082 const std::string feVolName;
9083 boost::shared_ptr<VolSideFe> sideOpFe;
9085 OpGetFieldGradientValuesOnSkinWithSingular(
9086 const std::string
field_name,
const std::string vol_fe_name,
9087 boost::shared_ptr<VolSideFe> side_fe)
9090 feVolName(vol_fe_name), sideOpFe(side_fe) {}
9095 if (
type != MBVERTEX)
9097 CHKERR loopSideVolumes(feVolName, *sideOpFe);
9102 boost::shared_ptr<VolSideFe> my_vol_side_fe_ptr =
9107 my_vol_side_fe_ptr->getOpPtrVector().push_back(
9109 "MESH_NODE_POSITIONS", data_hooke_element_at_pts->HMat));
9111 my_vol_side_fe_ptr->singularElement, my_vol_side_fe_ptr->singularDisp,
9113 data_hooke_element_at_pts->HMat));
9115 my_vol_side_fe_ptr->getOpPtrVector().push_back(
9117 my_vol_side_fe_ptr->detS,
9118 my_vol_side_fe_ptr->invSJac));
9120 my_vol_side_fe_ptr->getOpPtrVector().push_back(
9122 "SPATIAL_POSITION", data_hooke_element_at_pts->hMat));
9125 new OpGetFieldGradientValuesOnSkinWithSingular(
9126 "MESH_NODE_POSITIONS", fe_name.c_str(), my_vol_side_fe_ptr));
9132 new HookeElement::OpPostProcHookeElement<FaceElementForcesAndSourcesCore>(
9133 "MESH_NODE_POSITIONS", data_hooke_element_at_pts,
9153 post_proc->getOpPtrVector().push_back(
9158 post_proc->getOpPtrVector().push_back(
9161 post_proc->getOpPtrVector().push_back(
9167 boost::shared_ptr<moab::Interface> post_proc_mesh_ptr(
9169 boost::shared_ptr<std::vector<EntityHandle>> map_gauss_pts_ptr(
9171 post_proc->getOpPtrVector().push_back(
9173 post_proc_mesh_ptr, map_gauss_pts_ptr,
mwlsApprox));
9177 "mwlsApprox not allocated");
9182 approx_internal_stress) {
9186 post_proc->getOpPtrVector().push_back(
9191 post_proc->getOpPtrVector().push_back(
9194 post_proc->getOpPtrVector().push_back(
9200 boost::shared_ptr<moab::Interface> post_proc_mesh_ptr(
9202 boost::shared_ptr<std::vector<EntityHandle>> map_gauss_pts_ptr(
9204 post_proc->getOpPtrVector().push_back(
9206 post_proc_mesh_ptr, map_gauss_pts_ptr,
mwlsApprox,
9208 post_proc->commonData)));
9211 "mwlsApprox not allocated");
9218 ss <<
"out_skin_" << step <<
".h5m";
9225 ss <<
"out_spatial_" << step <<
".h5m";
9231 CHKERR post_proc->postProcMesh.tag_get_handle(
"SPATIAL_POSITION_GRAD",
9233 CHKERR post_proc->postProcMesh.tag_delete(
th);
9235 CHKERR post_proc->writeFile(ss.str());
9239 auto fe_post_proc_face_contact =
9240 boost::make_shared<PostProcFaceOnRefinedMesh>(
mField);
9241 CHKERR fe_post_proc_face_contact->generateReferenceElementMesh();
9243 CHKERR fe_post_proc_face_contact->addFieldValuesPostProc(
"LAMBDA_CONTACT");
9244 CHKERR fe_post_proc_face_contact->addFieldValuesPostProc(
9245 "SPATIAL_POSITION");
9246 CHKERR fe_post_proc_face_contact->addFieldValuesPostProc(
9247 "MESH_NODE_POSITIONS");
9250 fe_post_proc_face_contact);
9253 std::ostringstream stm;
9254 stm <<
"out_contact_surface_" << step <<
".h5m";
9256 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"out file %s\n",
9258 CHKERR fe_post_proc_face_contact->postProcMesh.write_file(
9271 std::ostringstream ostrm;
9272 ostrm <<
"out_contact_integ_pts_" << step <<
".h5m";
9274 CHKERR PetscPrintf(PETSC_COMM_WORLD,
"out file %s\n",
9277 "PARALLEL=WRITE_PART");
9290 if (dof_ptr->get()->getEntType() != MBVERTEX) {
9291 dof_ptr->get()->getFieldData() = 0;
9293 int dof_rank = dof_ptr->get()->getDofCoeffIdx();
9294 double &fval = dof_ptr->get()->getFieldData();
9297 fval = coords[dof_rank];
9309 MBVERTEX, dof_ptr)) {
9311 int dof_rank = dof_ptr->get()->getDofCoeffIdx();
9312 double fval = dof_ptr->get()->getFieldData();
9317 coords[dof_rank] = fval;
9318 if (dof_ptr->get()->getName() !=
field_name) {
9345 CHKERR moab.get_connectivity(*eit, conn, num_nodes,
false);
9347 CHKERR moab.get_coords(conn, num_nodes, coords);
9348 const double dir[3] = {coords[3] - coords[0], coords[4] - coords[1],
9349 coords[5] - coords[2]};
9350 double dof[3] = {0, 0, 0};
9352 for (
int dd = 0;
dd != 3;
dd++) {
9356 for (
int dd = 0;
dd != 3;
dd++) {
9362 const int idx = dit->get()->getEntDofIdx();
9364 dit->get()->getFieldData() = 0;
9367 cerr << **dit << endl;
9368 cerr << dof[idx] << endl;
9373 dit->get()->getFieldData() = dof[idx];
9396 mField,
"MESH_NODE_POSITIONS", *eit, dit)) {
9397 dit->get()->getFieldData() = 0;
9427 for (
int ll = 0; ll !=
nb_levels; ll++) {
9429 CHKERR moab.get_adjacencies(ents, 3,
false, adj_tets,
9430 moab::Interface::UNION);
9432 adj_tets = intersect(bit_tets, adj_tets);
9434 CHKERR moab.get_connectivity(adj_tets, ents,
true);
9437 CHKERR moab.get_adjacencies(adj_tets.subset_by_type(MBTET), 1,
false, ents,
9438 moab::Interface::UNION);
9439 CHKERR moab.get_adjacencies(adj_tets.subset_by_type(MBTET), 2,
false, ents,
9440 moab::Interface::UNION);
9441 ents.merge(adj_tets);
9447 bit,
BitRefLevel().set(), MBTET,
"bit2_at_crack_front.vtk",
"VTK",
"");
9449 bit,
BitRefLevel().set(), MBPRISM,
"bit2_at_crack_front_prisms.vtk",
9467 if (getInterface<CPMeshCut>()->getEdgesBlockSet())
9481 std::vector<int> surface_ids,
9482 const int verb,
const bool debug) {
9485 auto all_bits = []() {
return BitRefLevel().set(); };
9487 auto get_edges_block_set = [
this]() {
9488 return getInterface<CPMeshCut>()->getEdgesBlockSet();
9508 if (get_edges_block_set()) {
9526 const std::vector<int> surface_ids,
const int verb,
const bool debug) {
9537 ->writeEntitiesAllBitLevelsByType(
BitRefLevel().set(), MBTET,
9538 "all_bits.vtk",
"VTK",
"");
9558 for (
auto f : *fields_ptr) {
9559 if (
f->getName().compare(0, 6,
"LAMBDA") == 0 &&
9560 f->getName() !=
"LAMBDA_ARC_LENGTH" &&
9561 f->getName() !=
"LAMBDA_CONTACT" &&
9562 f->getName() !=
"LAMBDA_CLOSE_CRACK") {
9577 std::vector<std::string> fe_surf_proj_list) {
9592 fe_surf_proj_list,
false);
9594 dm_elastic, dm_material, bit1,
9597 "MATERIAL_FORCES",
QUIET);
9599 "SURFACE_PROJECTION", surface_ids,
9600 fe_surf_proj_list,
QUIET);
9602 "CRACK_SURFACE_AREA",
QUIET);
9634 all_ents = subtract(all_ents, all_ents.subset_by_type(MBENTITYSET));
9635 bit_ents = intersect(bit_ents, all_ents);
9636 all_ents = subtract(all_ents, bit_ents);
9639 for (Range::iterator eit = all_ents.begin(); eit != all_ents.end(); eit++) {
9640 for (RefEntity_multiIndex::iterator eiit =
9641 refined_ents_ptr->get<
Ent_mi_tag>().lower_bound(*eit);
9642 eiit != refined_ents_ptr->get<
Ent_mi_tag>().upper_bound(*eit);
9644 cerr << **eiit <<
" " << eiit->get()->getBitRefLevel() << endl;
9651 "", &out_meshset, 1);
9666 dofs_ptr->get<
Unique_mi_tag>().lower_bound(FieldEntity::getLoBitNumberUId(
9669 dofs_ptr->get<
Unique_mi_tag>().upper_bound(FieldEntity::getHiBitNumberUId(
9672 std::distance(dit, hi_dit) != 1)
9674 "Arc length lambda field not defined, or not unique");
9683 "Null pointer, probably field not found");
9685 const int rank =
fieldPtr->getNbOfCoeffs();
9689 "Is assumed that the number of field coefficients is smaller than 3"
9696 if (
rval != MB_SUCCESS) {
9698 tagName.c_str(), rank, MB_TYPE_DOUBLE,
tH, MB_TAG_CREAT | MB_TAG_SPARSE,
9699 &*def_val.data().begin());
9701 if (
F != PETSC_NULL) {
9709 if (
F != PETSC_NULL) {
9710 CHKERR VecRestoreArray(
F, &aRray);
9717 if (dofNumeredPtr->getEntType() != MBVERTEX)
9719 int local_idx = dofNumeredPtr->getPetscLocalDofIdx();
9720 if (local_idx == -1)
9723 int dof_rank = dofNumeredPtr->getDofCoeffIdx();
9724 int rank = dofNumeredPtr->getNbOfCoeffs();
9728 if (
F != PETSC_NULL) {
9729 tag_val[dof_rank] = aRray[local_idx];
9731 tag_val[dof_rank] = dofNumeredPtr->getFieldData();
9740 elementOrientation = +1;
9741 Range adj_side_elems;
9743 bit, &face, 1, 3, adj_side_elems, moab::Interface::INTERSECT,
QUIET);
9744 adj_side_elems = adj_side_elems.subset_by_type(MBTET);
9745 if (adj_side_elems.empty()) {
9746 Range adj_tets_on_surface;
9750 bit_tet_on_surface, &face, 1, 3, adj_tets_on_surface,
9751 moab::Interface::INTERSECT, 0);
9752 adj_side_elems.insert(*adj_tets_on_surface.begin());
9754 if (adj_side_elems.size() != 1) {
9755 adj_side_elems.clear();
9757 bit, &face, 1, 3, adj_side_elems, moab::Interface::INTERSECT,
9759 Range::iterator it = adj_side_elems.begin();
9760 for (; it != adj_side_elems.end(); it++) {
9762 CHKERR m_field.
get_moab().get_connectivity(&*it, 1, nodes,
true);
9763 PetscPrintf(PETSC_COMM_WORLD,
"Connectivity %lu %lu %lu %lu\n", nodes[0],
9764 nodes[1], nodes[2], nodes[3]);
9767 MPI_Comm_rank(m_field.
get_comm(), &rank);
9770 CHKERR m_field.
get_moab().create_meshset(MESHSET_SET, out_meshset);
9771 CHKERR m_field.
get_moab().add_entities(out_meshset, adj_side_elems);
9773 CHKERR m_field.
get_moab().write_file(
"debug_error.vtk",
"VTK",
"",
9778 "Expect 1 tet but is %u", adj_side_elems.size());
9781 if (side_elem != 0) {
9782 int side_number, sense, offset;
9783 CHKERR m_field.
get_moab().side_number(side_elem, face, side_number, sense,
9786 elementOrientation = -1;
9789 if (useProjectionFromCrackFront || contactFaces.contains(
Range(face, face))) {
9790 Tag th_interface_side;
9794 CHKERR m_field.
get_moab().tag_get_data(th_interface_side, &face, 1, &side);
9796 elementOrientation = -1;
9805 EntityHandle fe_ent = numeredEntFiniteElementPtr->getEnt();
9810 auto constrain_nodes = [&]() {
9812 std::vector<int> lambda_dofs(9, -1);
9813 std::vector<double> lambda_vals(9, 0);
9818 auto row_dofs = getRowDofsPtr();
9819 auto lo_uid_lambda = DofEntity::getLoFieldEntityUId(
9820 getFieldBitNumber(lambdaFieldName), get_id_for_min_type<MBVERTEX>());
9821 auto hi_uid_lambda = DofEntity::getHiFieldEntityUId(
9822 getFieldBitNumber(lambdaFieldName), get_id_for_max_type<MBVERTEX>());
9824 for (
auto it = row_dofs->lower_bound(lo_uid_lambda);
9825 it != row_dofs->upper_bound(hi_uid_lambda); ++it) {
9826 int side_number = it->get()->getSideNumberPtr()->side_number;
9827 if (side_number > 2)
9829 const int dof_idx = 3 * side_number + it->get()->getEntDofIdx();
9830 lambda_dofs[dof_idx] = it->get()->getPetscGlobalDofIdx();
9831 lambda_vals[dof_idx] = it->get()->getFieldData();
9836 "Only nine Lagrange multipliers can be on element (%d)", sum_dd);
9838 std::vector<int> positions_dofs(18, -1);
9839 std::vector<double> positions_vals(18, 0);
9840 std::vector<bool> positions_flags(18,
true);
9842 auto lo_uid_pos = DofEntity::getLoFieldEntityUId(
9843 getFieldBitNumber(spatialFieldName), get_id_for_min_type<MBVERTEX>());
9844 auto hi_uid_pos = DofEntity::getHiFieldEntityUId(
9845 getFieldBitNumber(spatialFieldName), get_id_for_max_type<MBVERTEX>());
9848 for (
auto it = row_dofs->lower_bound(lo_uid_pos);
9849 it != row_dofs->upper_bound(hi_uid_pos); ++it) {
9850 int side_number = it->get()->getSideNumberPtr()->side_number;
9851 int dof_idx = 3 * side_number + it->get()->getEntDofIdx();
9852 positions_dofs[dof_idx] = it->get()->getPetscGlobalDofIdx();
9853 positions_vals[dof_idx] = it->get()->getFieldData();
9854 Range vert =
Range(it->get()->getEnt(), it->get()->getEnt());
9855 if (masterNodes.contains(vert)) {
9856 positions_flags[dof_idx] =
false;
9860 const double beta = area * aLpha;
9862 case CTX_SNESSETFUNCTION:
9863 CHKERR VecSetOption(snes_f, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
9864 for (
int ii = 0; ii != 9; ++ii) {
9865 if (lambda_dofs[ii] == -1)
9867 double gap = positions_vals[0 + ii] - positions_vals[9 + ii];
9868 double val1 = beta * gap;
9869 CHKERR VecSetValue(snes_f, lambda_dofs[ii], val1, ADD_VALUES);
9870 double val2 = beta * lambda_vals[ii];
9871 if (positions_flags[0 + ii]) {
9872 CHKERR VecSetValue(snes_f, positions_dofs[0 + ii], +val2, ADD_VALUES);
9874 if (positions_flags[9 + ii]) {
9875 CHKERR VecSetValue(snes_f, positions_dofs[9 + ii], -val2, ADD_VALUES);
9879 case CTX_SNESSETJACOBIAN:
9880 for (
int ii = 0; ii != 9; ++ii) {
9881 if (lambda_dofs[ii] == -1)
9883 CHKERR MatSetValue(snes_B, lambda_dofs[ii], positions_dofs[0 + ii],
9884 +1 * beta, ADD_VALUES);
9885 CHKERR MatSetValue(snes_B, lambda_dofs[ii], positions_dofs[9 + ii],
9886 -1 * beta, ADD_VALUES);
9887 if (positions_flags[0 + ii]) {
9888 CHKERR MatSetValue(snes_B, positions_dofs[0 + ii], lambda_dofs[ii],
9889 +1 * beta, ADD_VALUES);
9891 if (positions_flags[9 + ii]) {
9892 CHKERR MatSetValue(snes_B, positions_dofs[9 + ii], lambda_dofs[ii],
9893 -1 * beta, ADD_VALUES);
9903 auto constrain_edges = [&]() {
9906 map<EntityHandle, EntityHandle> side_map;
9909 for (
auto s : {0, 1, 2}) {
9928 map<EntityHandle, std::vector<int>> side_map_lambda_dofs;
9929 map<EntityHandle, std::vector<double>> side_map_lambda_vals;
9930 map<EntityHandle, std::vector<int>> side_map_positions_dofs;
9931 map<EntityHandle, std::vector<double>> side_map_positions_vals;
9932 map<EntityHandle, double> side_map_positions_sign;
9934 for (
auto e : ents) {
9936 auto row_dofs = getRowDofsPtr();
9937 auto lo_uid_lambda = row_dofs->lower_bound(DofEntity::getLoFieldEntityUId(
9938 getFieldBitNumber(lambdaFieldName), e));
9939 auto hi_uid_lambda = row_dofs->upper_bound(DofEntity::getHiFieldEntityUId(
9940 getFieldBitNumber(lambdaFieldName), e));
9941 const int nb_lambda_dofs = std::distance(lo_uid_lambda, hi_uid_lambda);
9943 std::vector<int> lambda_dofs(nb_lambda_dofs, -1);
9944 std::vector<double> lambda_vals(nb_lambda_dofs, 0);
9946 for (
auto it = lo_uid_lambda; it != hi_uid_lambda; ++it) {
9947 const int dof_idx = it->get()->getEntDofIdx();
9948 lambda_dofs[dof_idx] = it->get()->getPetscGlobalDofIdx();
9949 lambda_vals[dof_idx] = it->get()->getFieldData();
9952 if (nb_lambda_dofs) {
9953 side_map_lambda_dofs[e] = lambda_dofs;
9954 side_map_lambda_vals[e] = lambda_vals;
9955 side_map_lambda_dofs[side_map[e]] = lambda_dofs;
9956 side_map_lambda_vals[side_map[e]] = lambda_vals;
9959 auto lo_uid_pos = row_dofs->lower_bound(DofEntity::getLoFieldEntityUId(
9960 getFieldBitNumber(spatialFieldName), e));
9961 auto hi_uid_pos = row_dofs->upper_bound(DofEntity::getHiFieldEntityUId(
9962 getFieldBitNumber(spatialFieldName), e));
9963 const int nb_dofs = std::distance(lo_uid_pos, hi_uid_pos);
9965 std::vector<int> positions_dofs(nb_dofs, -1);
9966 std::vector<double> positions_vals(nb_dofs, 0);
9967 side_map_positions_sign[e] = (!nb_lambda_dofs) ? 1 : -1;
9970 for (
auto it = lo_uid_pos; it != hi_uid_pos; ++it) {
9971 int dof_idx = it->get()->getEntDofIdx();
9972 positions_dofs[dof_idx] = it->get()->getPetscGlobalDofIdx();
9973 positions_vals[dof_idx] = it->get()->getFieldData();
9976 side_map_positions_dofs[e] = positions_dofs;
9977 side_map_positions_vals[e] = positions_vals;
9980 for (
auto m : side_map_positions_dofs) {
9983 if (side_map_lambda_dofs.find(e) != side_map_lambda_dofs.end()) {
9985 auto &lambda_dofs = side_map_lambda_dofs.at(e);
9986 auto &lambda_vals = side_map_lambda_vals.at(e);
9987 auto &positions_dofs = side_map_positions_dofs.at(e);
9988 auto &positions_vals = side_map_positions_vals.at(e);
9989 auto sign = side_map_positions_sign.at(e);
9991 const double beta = area * aLpha;
9994 case CTX_SNESSETFUNCTION:
9995 for (
int ii = 0; ii != lambda_dofs.size(); ++ii) {
9996 double val1 = beta * positions_vals[ii] *
sign;
9997 CHKERR VecSetValue(snes_f, lambda_dofs[ii], val1, ADD_VALUES);
9999 for (
int ii = 0; ii != positions_dofs.size(); ++ii) {
10000 double val2 =
sign * beta * lambda_vals[ii];
10001 CHKERR VecSetValue(snes_f, positions_dofs[ii], +val2, ADD_VALUES);
10004 case CTX_SNESSETJACOBIAN:
10005 for (
int ii = 0; ii != lambda_dofs.size(); ++ii) {
10006 CHKERR MatSetValue(snes_B, lambda_dofs[ii], positions_dofs[ii],
10007 beta *
sign, ADD_VALUES);
10009 for (
int ii = 0; ii != positions_dofs.size(); ++ii) {
10010 CHKERR MatSetValue(snes_B, positions_dofs[ii], lambda_dofs[ii],
10011 beta *
sign, ADD_VALUES);
10023 CHKERR constrain_nodes();
10024 CHKERR constrain_edges();
10036 constexpr
double alpha = 1.e-5;
10041 double temp = 250.;
10042 double z = t_coords(2);
10043 if ((-10. < z && z < -1.) || std::abs(z + 1.) < 1e-15) {
10044 temp = 10. / 3. * (35. - 4. * z);
10046 if ((-1. < z && z < 2.) || std::abs(z - 2.) < 1e-15) {
10047 temp = 10. / 3. * (34. - 5. * z);
10049 if (2. < z && z < 10.) {
10050 temp = 5. / 4. * (30. + 17. * z);
10053 t_thermal_strain(
i,
j) = alpha * (
temp - 250.) *
t_kd(
i,
j);
10054 return t_thermal_strain;