v0.14.0
Loading...
Searching...
No Matches
VolumeElementForcesAndSourcesCoreOnContactPrismSide.cpp
Go to the documentation of this file.
1/** \file VolumeElementForcesAndSourcesCoreOnContactPrismSide.cpp
2
3\brief Implementation of volume element on side
4
5*/
6
7
8
9namespace MoFEM {
10
14
18
19 if (!sidePtrFE)
20 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
21 "Side volume element not set");
22
23 const EntityHandle volume_entity =
25
26 Range adj_faces;
27 CHKERR sidePtrFE->mField.getInterface<BitRefManager>()->getAdjacenciesAny(
28 volume_entity, 2, adj_faces);
29 SideNumber_multiIndex &side_table = const_cast<SideNumber_multiIndex &>(
30 numeredEntFiniteElementPtr->getSideNumberTable());
31
32 SideNumber_multiIndex &side_volume_table =
33 const_cast<SideNumber_multiIndex &>(
34 sidePtrFE->numeredEntFiniteElementPtr->getSideNumberTable());
35
36 adj_faces = adj_faces.subset_by_type(MBTRI);
37 bool face_common = false;
38 SideNumber_multiIndex::nth_index<0>::type::iterator sit;
39 SideNumber_multiIndex::nth_index<0>::type::iterator sit_volume;
40 EntityHandle face_ent;
41 for (Range::iterator it_tris = adj_faces.begin(); it_tris != adj_faces.end();
42 it_tris++) {
43 sit = side_table.get<0>().find(*it_tris);
44 sit_volume = side_volume_table.get<0>().find(*it_tris);
45 if (sit != side_table.get<0>().end()) {
46 face_ent = *it_tris;
47 face_common = true;
48 break;
49 }
50 }
51 if (!face_common) {
52 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Error no common tets");
53 }
54
55 auto contact_prism_ptr_fe =
57
58 faceSense = (*sit)->sense;
59 faceSideNumber = (*sit)->side_number;
60 int side_of_vol_number = (*sit_volume)->side_number;
61 const EntityHandle *ent_on_vol;
62 int num_nodes = 3;
63 CHKERR contact_prism_ptr_fe->mField.get_moab().get_connectivity(
64 face_ent, ent_on_vol, num_nodes, true);
65
66 fill(tetConnMap.begin(), tetConnMap.end(), -1);
67 for (int nn = 0; nn != 3; ++nn) {
68 faceConnMap[nn] = std::distance(conn, find(conn, &conn[4], ent_on_vol[nn]));
69
70 tetConnMap[faceConnMap[nn]] = nn;
71 if (faceConnMap[nn] > 3) {
72 cerr << " faceConnMap[nn] " << faceConnMap[nn] << "\n";
73 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
74 "No common node on face and element can not be found");
75 }
76 }
77
78 oppositeNode = std::distance(tetConnMap.begin(),
79 find(tetConnMap.begin(), tetConnMap.end(), -1));
80
81 const int nb_gauss_pts =
82 contact_prism_ptr_fe->getGaussPtsMasterFromEleSide().size2();
83 gaussPts.resize(4, nb_gauss_pts, false);
84 gaussPts.clear();
85
86 boost::shared_ptr<EntitiesFieldData> dataH1_on_face;
87
88 if (side_of_vol_number == 3) {
89 dataH1_on_face = boost::shared_ptr<EntitiesFieldData>(
90 contact_prism_ptr_fe->getDataOnMasterFromEleSide()[H1]);
91 } else if (side_of_vol_number == 4) {
92 dataH1_on_face = boost::shared_ptr<EntitiesFieldData>(
93 contact_prism_ptr_fe->getDataOnSlaveFromEleSide()[H1]);
94 } else {
95 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
96 "Side element for contact: Wrong face for contact!");
97 }
98
99 const MatrixDouble &face_shape_funtions =
100 dataH1_on_face->dataOnEntities[MBVERTEX][0].getN(NOBASE);
101
102 const double tet_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
103 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
104 gaussPts(0, gg) =
105 face_shape_funtions(gg, 0) * tet_coords[3 * faceConnMap[0] + 0] +
106 face_shape_funtions(gg, 1) * tet_coords[3 * faceConnMap[1] + 0] +
107 face_shape_funtions(gg, 2) * tet_coords[3 * faceConnMap[2] + 0];
108 gaussPts(1, gg) =
109 face_shape_funtions(gg, 0) * tet_coords[3 * faceConnMap[0] + 1] +
110 face_shape_funtions(gg, 1) * tet_coords[3 * faceConnMap[1] + 1] +
111 face_shape_funtions(gg, 2) * tet_coords[3 * faceConnMap[2] + 1];
112 gaussPts(2, gg) =
113 face_shape_funtions(gg, 0) * tet_coords[3 * faceConnMap[0] + 2] +
114 face_shape_funtions(gg, 1) * tet_coords[3 * faceConnMap[1] + 2] +
115 face_shape_funtions(gg, 2) * tet_coords[3 * faceConnMap[2] + 2];
116 if (side_of_vol_number == 3) {
117 gaussPts(3, gg) =
118 contact_prism_ptr_fe->getGaussPtsMasterFromEleSide()(2, gg);
119 } else {
120 gaussPts(3, gg) =
121 contact_prism_ptr_fe->getGaussPtsSlaveFromEleSide()(2, gg);
122 }
123 }
125}
126
127} // namespace MoFEM
@ NOBASE
Definition definitions.h:59
@ H1
continuous field
Definition definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
constexpr int order
multi_index_container< boost::shared_ptr< SideNumber >, indexed_by< ordered_unique< member< SideNumber, EntityHandle, &SideNumber::ent > >, ordered_non_unique< composite_key< SideNumber, const_mem_fun< SideNumber, EntityType, &SideNumber::getEntType >, member< SideNumber, signed char, &SideNumber::side_number > > > > > SideNumber_multiIndex
SideNumber_multiIndex for SideNumber.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
Managing BitRefLevels.
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
ForcesAndSourcesCore * sidePtrFE
Element to integrate on the sides.
MatrixDouble gaussPts
Matrix of integration points.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.